在多重檢定，例如變異數分析 (ANOVA) 中的多重比較，一直被關注的是族錯誤率 (FWER, family-wise error rate)，但在控制族錯誤率的同時，不免造成檢定力 (power of a test) 的低下。為此，Benjamini and Hochberg (1995) 提出了「偽發現率 (FDR, False Discovery Rate)」的概念及相應的多重檢定程序。
偽發現率，定義為
在所有被拒絕的虛無假說中，其實犯了型Ⅰ誤者所佔比例的期望值。
如果檢定問題族 Hi vs. Ki, i = 1, ..., m, 之中有 r 個虛無假說 Hi 被拒絕，但其中有 k 個 Hi 其實是成立的，則犯了型Ⅰ誤者佔拒絕之 Hi 比例是 k/r。不過，此處 k 和 r 都是隨機量，因此 FDR 以其期望值定義。
以疾病篩檢來看，假設敏感敏感度是 β，相當於假說檢定中的檢定力；特異度是 q，相當於型Ⅰ誤機率 α 的 1-補數，即 q = 1-α；犯病率是 d。則 FDR 接近
期望偽陽性數 / 期望總陽性數
= (1-d)(1-q)/[(1-d)(1-q)+dβ] = (1-d)α/((1-d)α + dβ]
上式是先求平均再算比例，FDR 則是計算比例的平均，因此只說「接近」。事實上 FDR 應略高於上式，但因拒絕數可能為 0，FDR 的正式定義可能使「略高」這結論不成立。前列算式顯示在疾病篩檢中，檢驗結果是陽性者包括真陽性（有病且檢驗結果是陽性）與偽陽性（無病但檢驗結果是陽性）。在統計假說檢定，被拒絕的 Hi 也分二類：不該被拒絕的 (Hi 成立），與應被拒絕的 (Ki 成立）。與疾病篩檢不同的是：在疾病篩檢問題，敏感度與特異度可以假設是常數；在假說檢定問題，Hi 不成立者其檢定力固然不是常數；Hi 成立者其型Ⅰ誤機率也不一定是常數，但如果把問題簡化可以其上限，顯著水準代之，如果各檢定問題，Hi vs. Ki，都採相同固定水準檢定，勉強可當作是常數。假設 m = 1, 只有一個檢定，H 對 K，當 H 成立時 (d = 0)，依前列計算式得 1，當 H 不成立時得值 0。
假設 m 個 Hi 中，成立的有 n 個，所以 m-n 個 Ki 成立。在 n 個成立的 Hi 中，U 個不被拒絕，V 個被拒絕；而 m-n 個不成立的 Hi 中，T 個不被拒絕，S 個被拒絕。故拒絕數是 V + S，而 FDR 定義為
FDR = E[V/(V+S)] = E[V/(V+S) | V+S>0] P[V+S > 0]
即 V+S = 0 時 V/(V+S) 被定義為 0。設 m = 1, 僅有一個檢定問題，H 成立時 S = T = 0, FDR 等於 P[V+S>0] = P[V = 1] 是型Ⅰ誤機率；K 成立時 U = V = 0, FDR 等於 0。假設 m 個檢定問題的檢定統計量是相互獨立的，又假設 m 個檢定都是 α 水準檢定，為簡單計更假設型Ⅰ誤機率都是 α，於是 V 是 Bin(n, α) 隨機變數，而 S 則是一些 Bernoulli(β_i) 變數的和，則
P[V + S > 0] = 1 - (1-α)^n Π" (1-β_i)
E[V] = nα, E[S] = Σ" β_i
以上 Π" 和 Σ" 分別是在 Ki 成立之範圍內的 i 做連乘和連加，p_i 是 K_i 成立之檢定的檢定力，所以
FDR ≒ E[V]/E[V+S] = nα/(nα + Σ" β_i)
當 (1-α)^n Π" (1-β_i) 很小時，右式其實略小於 FDR, 如前面所述。由上式看來：若所有 Hi 都成立，則 FDR 接近 １。不過，如果 P[V > 0] 並不接近 1，如同 m = 1 的情形，上列 FDR 「近似式」並不準確，或者說，並不適用。事實上所有 Hi 都成立時，所有的拒絕 Hi 都是型Ⅰ錯誤，也就是所有陽性都是偽陽性，因此
FDR = 1．P[V + S > 0] = P[V > 0] = FWER
在前面相互獨立及同為 α 水準檢定的假設下，其值是 1 - (1-α)^m ≒ 1 - e^{-mα}。但上列 FDR = FWER 是定義上的等值，不論這些檢定是否相互獨立，是否同一個顯著水準。
由「近似式」來看，FDR 的大小很大比例取決於所考慮的檢定問題族有多少 Ki 是成立的，以及對應的檢定力大小，例如另一個極端情形，所有 Hi 都不成立，也就是說都是 Ki 成立，則 FDR 是 0。正確的 FDR 即使在前述相互獨立及同為 α 水準檢定的假設下也不容易計算，而「近似式」可能不準確，但可能可以設法改進。令 f(V,S) = V/(V+S), 則 f(V,S) 在 (μ_v, μ_s) 做 Taylor 展開得
f(V,S) = μ_v/(μ_v+μ_s)
+ (V-μ_v)μ_s/(μ_v+μ_s)^2 - (S-μ_s)μ_v/(μ_v+μ_s)^2
- (V-μ_v)^2 μ_s/(μ_v+μ_s)^3 + (S-μ_s)^2μ_v/(μ_v+μ_s)^3
+ (1/2) (V-μ_v)(S-μ_s)(μ_v-μ_s)/(μ_v+μ_s)^3
+ (V-μ_v)^3 μ_s/(μ_v+μ_s)^4 - (S-μ_s)^3μ_v/(μ_v+μ_s)^4
- (1/3) (V-μ_v)^2(S-μ_s)(μ_v-2μ_s)/(μ_v+μ_s)^4
- (1/3) (V-μ_v)(S-μ_s)^2(2μ_v-μ_s)/(μ_v+μ_s)^4
+ ...
在前述相互獨立的及同為 α 水準檢定的條件下，取 2 階近似得
FDR ≒ μ_v/(μ_v+μ_s) - Var(V)μ_s/(μ_v+μ_s)^3 + Var(S)μ_v/(μ_v+μ_s)^3
= nα/(nα + Σ" β_i) - nα(1-α)(Σ" β_i)/(nα + Σ" β_i)^3
+ Σ" β_i(1-β_i) nα/(nα + Σ" β_i)^3
第二階項是 O(1/m) 等級，也就是說如果 m 較大，改進項效果較小，同時也暗示更高階項的不重要。
Benjamini and Hochberg (1995) 提出一個可以控制 FDR 的檢定程序（依該論文的備註 (remark)，此程序其實 Simes (1988) 已提到）：
0. 設定一個標準，記為 q*。
1. 將各檢定問題的 p-值排序,假設 p_1 < ... < p_m。
2. 令 k = max{i: p_i ≦ (i/m)q*}, 拒絕虛無假說 Hi, i=1,...,k。
由於 p_i > (i/m)q* 並不保證 j > i 時 p_j > (j/m)q*，因此上列程序所拒絕的 Hi 可能包含不滿足其 p_i ≦ (i/m)q* 者，也就是說：依上列程序，被拒絕的 Hi （依其 p-值排序）包括：
(a) p_i ≦ (i/m)q* 所對應的 Hi；
(b) 存在較大的 p_j 滿足 p_j ≦ (j/m)q* 所對應的 Hi。
如果各檢定的檢定統計量都是連續的，如果所有虛無假說都是成立的，所有 p-值都服從 0-1 之間的均勻分布，排序過的 p_i 等於原來 m 個 p-值的順序統計量，為了方便，改用 U＿i 取代排序後的 p-值，故
1 - FDR = 1 - FWER = P{ k = 0 }
= P[ U_i > (i/m)q*, i = 1, ..., m ]
= ∫_[q*,1]…∫_[(i/m)q*,u_{i+1}]…∫_[q*/m,u_2] m! du_1…du_i…du_m
= 1 - q*
因此，Benjamini and Hochberg 程序控制了 FWER 不超過 q*。不過所謂FWER 不超過 q* 只是指所有 m 個虛無假說都成立的情況，如果 m 個 Hi 中有些是不成立的，那麼其餘成立的 Hi 會有至少一個被錯誤的拒絕的機率仍能控制在 q* 之內嗎？Benjamini and Hochberg 引述 Hommel (1088) 的結果是：不能。例如，假設 n = m-1, q* = α，諸 Hi, i = 1, ..., n 是成立的，唯 Hm 不成立。仿「Holm 的 Bonferroni 校正」一文末段的計算，設 p* = min{p_1, ..., p_m}，p° = p_m，在 m 個檢定統計量相互獨立且均是連續型的條件下，給定 p°, 諸虛無假說 Hi 中至少一個被錯誤地拒絕的機率，當 p° ≦ α/m 時是
P[p* < p° | p°] + P[p* ≧ p° 且 p* ≦ 2α/m | p°] = 1 - (1-2α/m)^{m-1} ≒ [2(m-1)/m]α
當 p° > α/m 時得
P[p* < p° 且 p* ≦ α/m | p°] = 1 - [1 - α/m]^{m-1} ≒ [(m-1)/m]α
故
FWER ≒ P[p° ≦ α/m] [2(m-1)/m]α + P[p° > α/m] [(m-1)/m]α
≧ [(m-1)/m]α + [(m-1)/m^2]α^2
雖然最後一式比 α 略小，但它只是 FWER 的下界，而真正的 FWER 把更大權重放在 [2(m-1)/m]α，只要 m > 2 它就比 α 來得大。因此，Benjamini and Hochberg 程序並不能在任意情形嚴格控制 FWER。話說回來，此程序意在控制 FDR，它不能嚴格意義上控制 FWER 就不那麼重要了。
由於 FDR = E[V/(V+S)] = E[V/(V+S) | V+S>0] P[V+S > 0]，當諸虛無假說 Hi 中不全成立時，一是 Hi 不成立對應的 p-值傾向於較小，Hi 被拒絕的機率較大；二是 Hi 不成立的多則 S 大 V 小，因此 V/(V+S) 也小。不過有 Hi 不成立使 P[V+S > 0] 增大，甚至可能大很多，因此從 FDR 的定義式倒是難以定論在有 Hi 不成立時 FDR 是否比全部 Hi 都成立時來得小。如前，令 m 個 Hi 中有 n 個是成立的，當 n = m 也就是所有 Hi 都成立，則如前述，在適當條件 (各檢定統計量均為連續型且相互獨立）下，Benjamini and Hochberg 程序（簡稱 BH 程序）控制了 FWER 不超過 q*。當 n = 0，所有被拒絕的 Hi 都是不成立的虛無假說，所以 FDR = 0。在 Benjamini and Hochberg (1995) 文中有個引理：
[引理] 若 m 個檢定的統計量相互獨立，在給定對應不成立的 Hi 的 p-值是 p'(1), ..., p'(n')，其中 n' = m - n，條件下。
E[ V/(V+S) | p'(1), ..., p'(n')] ≦ (n/m) q*
既然 V/(V+S) 的條件期望值都不超過 (n/m)q*, 其無條件期望值當然也不超過 (n/m) q*。這也就是說：成立的虛無假說佔比愈大，則 FDR 愈大；一個附帶的結論是 BH 程序確實把 FDR 控制在 q* 以下。
當 m = 1 也就是只有一個檢定 H 對 K 時，H 成立則 FDR = q*，K 成立則 FDR = 0，引理結論自然成立。假設上列引理結論在 m ≦ m° 時都成立，如果能證明當 m = m° + 1 時引理結論也成立，依數學歸納法原理，就證明了引理成立。前面說了，n = 0 時 FDR = 0，所以只需看 n > 0 的情形，事實上只需看 0 < n < m 的情形，因為 n = m 即所有 Hi 都成立，FDR = FWER ≦ q*。
假設對應不成立的 Hi 的 p-值是 p'(1)<…< p'(n')；對應成立的 Hi 的 p-值不排序是 p_1, …, p_n，又令 p* = max{p_1, …, p_n}。注意
E[ V/(V+S) | p'(1), ..., p'(n')]
= E[E[ V/(V+S) | p*, p'(1), ..., p'(n')] | p'(1), ..., p'(n')]
令 j° = max{j: p'(j) ≦ [(n+j)/m] q*}, p° = [(n+j°)/m]q*。假設檢定統計量是連績型，則
E[ V/(V+S) | p'(1), ..., p'(n')]
= ∫_[0,p°] [E[ V/(V+S) | p*, p'(1), ..., p'(n')] n p^{n-1} dp
+ ∫_[p°,1] [E[ V/(V+S) | p*, p'(1), ..., p'(n')] n p^{n-1} dp
式中 np^{n-1} 是 p* 的機率密度。當 p* ≦ p° 時，p'(j°) 是所有 m 個 p-值中第 n+j° 個，因此依 BH 程序，n 個成立的虛無假說及 j° 個不成立的虛無假說都被拒絕，
E[ V/(V+S) | p*, p'(1), ..., p'(n')] = n/(n+j°)
故
∫_[0,p°] [E[ V/(V+S) | p*, p'(1), ..., p'(n')] n p^{n-1} dp
= [n/(n+j°)] p°^n = (n/m) q* p°^{n-1}
當 p* > p° 時，如果,p'(j°) ≦ p° < p* < p'(j°+1)，則 p* 是所有 m 個 p-值中第 n+j° 個，而它大於 p° = [(n+j°)/m] q*, 因此對應的虛無假說不被拒絕。而若存在 j > j° 使
p'(j°) ≦ p° ＜ p'(j) ≦p* < p'(j+1)
則 p* 是 m 個 p-值中第 n+j 個，p* > p'(j) > [(n+j)/m] q*，故 p* 對應的虛無假說也不被拒絕。但 p* 對應的虛無假說，及之後（前面的 p'(j°+1), p'(j°+2) 等及 p'(j+1), p'(j+2) 等）對應的假說不被拒絕，則
E[ V/(V+S) | p*, p'(1), ..., p'(n')]
可看成 m-j-1 個 Hi 其中 n-1 個成立，j 個不成立，對應 p-值分別是 p_i/p*, p'(j)/p*, BH 程序用
[k/(n+j-1)] [(n+j-1)/m] (q*/p*)
做比較時的
E[ V/(V+S) | p'(1)/p*, ..., p'(n')/p*]
依假設，在 m = n+j-1 ≦ m° 時上列條件期望值不超過
[(n-1)/(n+j-1)] [(n+j-1)/m] (q*/p*) = [(n-1)/m] (q*/p*)
故
∫_[p°,1] [E[ V/(V+S) | p*, p'(1), ..., p'(n')] n p^{n-1} dp
≦ ∫_[p°,1] [(n-1)/m](q*/p) n p^{n-1} dp
= [n/(n-1)] [(n-1)/m] q* (1-p°^{n-1} = (n/m)q* (1-p°^{n-1}
將 E[ V/(V+S) | p'(1), ..., p'(n')] 兩段積分相加，結果得此條件期望值上限 (n/m) q* 在 m = m° + 1 時成立，引理得證。
前述「引理」的前題是所有 m 個檢定相互獨立，雖然實際上 Benjamini and Hochberg 特別強調那些 Ki 成立對應的檢定之間不必要求相互獨立，但我們並不知 m 個檢定問題中有幾個 Hi 是成立的，更遑論明確知道哪些個 Hi 成立，哪些個 Ki 成立，因此只要求諸 Hi 成立對應的檢定之間相互獨立是不合宜的。同時這也指出一個問題：如果檢定族 Hi vs. Ki, i = 1, ..., n 並非相互獨立，前面的證明並不適用。如此，「引理」的結論是否仍然成立？如前，令 p* = max{p_1, …, p_n}，而 p° = [(n+j°)/m]q* 其中 j° = max{j: p'(j) ≦ [(n+j)/m] q*}。令 F(p) 為 p* 在給定 p'(1), ..., p'(n') 值之後的條件分布函數，則
E[ V/(V+S) | p'(1), ..., p'(n')]
= ∫_[0,p°] [E[ V/(V+S) | p*, p'(1), ..., p'(n')] dF(p*)
+ ∫_[p°,1] [E[ V/(V+S) | p*, p'(1), ..., p'(n')] dF(p*)
= n/(n+j°) F(p°) + ∫_[p°,1] [E[ V/(V+S) | p*, p'(1), ..., p'(n')] dF(p*)
≦ n/m F(p°)/p° + ∫_[p°,1] [(n-1)/m](q*/p*) dF(p*)
= (n/m)q*{F(p°)/p° + ∫_[p°,1] [(n-1)/n](1/p*) dF(p*)}
因此，在一般情形「引理」的結論是否能延沿用先前的數學歸納證法得證，關鍵在於是否
F(p°)/p° + ∫_[p°,1] [(n-1)/n](1/p*) dF(p*) ≦ 1
留言列表
統計 (9)
{{ article.title }}