3つ以上のグループの全ての組み合わせを比較するとき、Tukey法、Bonferroni法、Benjamini-Hochberg法 (BH法) の結果がどれくらい違うのか数値計算で調べてみました。
結論から書くと、調べた条件下で陽性判定の出やすさは、補正なし > BH法 > Tukey法 & Bonferroni法でした。Tukey法とBonferroni法はほとんど同じでした。
比較条件
| 項目 | 値 |
|---|---|
| グループ数 | 3, 4, 5, 6 |
| 各グループのデータ点の数 | 全て100 |
| 分布 | 正規分布 |
| 分散 | 全て1 |
| 平均 | 各グループに0からdまでを等間隔に割り当てる、dは0から1まで |
| サンプリング法 | ランダムでなく、分位関数から等間隔に選択 |
| 有意水準 or FDR | 0.05 |
結果
上から順に、補正なし、BH法、Tukey法、Bonferroni法です。横軸はグループ間の平均の差の最大値、縦軸はp値 (正確にはBH法ではq値) です。赤線は0.05の水準を表します。
各グループの平均は等間隔に設定したため、線が重なっていることに注意して下さい。差の大きい順 (p値の小さい順) に以下のようにまとまっています。
- 3グループ
- (1,3)
- (1,2), (2,3)
- 4グループ
- (1,4)
- (1,3), (2,4)
- (1,2), (2,3), (3,4)
- 5グループ
- (1,5)
- (1,4), (2,5)
- (1,3), (2,4), (3,5)
- (1,2), (2,3), (3,4), (4,5)
- 6グループ
- (1,6)
- (1,5), (2,6)
- (1,4), (2,5), (3,6)
- (1,3), (2,4), (3,5), (4,6)
- (1,2), (2,3), (3,4), (4,5), (5,6)
赤線を横切る位置 (どれだけ平均差があれば有意となるか) に着目すると、陽性判定の出やすさが補正なし > BH法 > Tukey法 & Bonferroni法の順になっていることが分かります。BH法では2番目以降のものがTukey法やBonferroni法よりも早く赤線を下回っています。
ちなみに、BH法では、p値が最大のものは補正なしと同じ、p値が最小のものはBonferroni法と同じです。p値の順位によって掛け算される値が違うため、d=1から左に見ていくと途中で衝突します。
考察
他の条件ではまた違うかもしれませんが、今回の結果をみる限りはTukey法の代わりにBonferroni法を使っても実用上あまり問題ないかもしれません。少なくとも、Tukey法は正規性と等分散性を仮定しているため、データによっては使ってよいかどうか判断に迷うことがありますが、Bonferroni法なら気にせず使えて便利という面はあります1。
もちろん、この結果からBH法が良いと考えることも可能です。ただし、FWER制御とFDR制御の違いに注意して下さい。つまり、陽性と判定されたもののうち5%程度は誤検出かもしれないがそれでも構わない、という場合にBH法が使えます。
ソースコード
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
sns.set_context('talk', font_scale=0.8)
def calculate_q(p_seq):
p_arr = np.asarray(p_seq)
N = len(p_arr)
idx_arr = np.argsort(p_arr)
q_arr = p_arr[idx_arr] * N / (np.arange(N) + 1)
q_arr = np.minimum.accumulate(q_arr[::-1])[::-1]
q_arr[idx_arr] = q_arr.copy()
return q_arr
def test_multiple_comparision2():
n = 100
n_group_list = [3, 4, 5, 6]
d_arr = np.linspace(0, 1, 50)
method_list = ['raw', 'bh', 'tukey', 'bonferroni']
# NOTE: ppf(0)=-inf and ppf(1)=inf are removed
data = stats.norm().ppf(np.linspace(0,1,n+2)[1:-1])
def calc_tukey(Y, triu):
return stats.tukey_hsd(*Y.T).pvalue[triu]
def calc_raw(Y, triu):
return stats.ttest_ind(Y[:,triu[0]], Y[:,triu[1]])[1]
# loop
out_list = []
for n_group in n_group_list:
X = np.tile(data.reshape(-1,1), n_group)
triu = np.triu_indices(n_group, 1)
for d in d_arr:
Y = X + np.linspace(0, d, n_group)
for method in method_list:
if method == 'tukey':
p_arr = calc_tukey(Y, triu)
elif method == 'raw':
p_arr = calc_raw(Y, triu)
elif method == 'bonferroni':
p_arr = np.clip(calc_raw(Y, triu) * triu[0].size, 0, 1)
elif method == 'bh':
p_arr = calculate_q(calc_raw(Y, triu))
df_tmp = pd.DataFrame()
df_tmp['p_val'] = p_arr
df_tmp['method'] = method
df_tmp['n_group'] = n_group
df_tmp['d'] = d
out_list.append(df_tmp)
df = pd.concat(out_list, axis=0)
# plot
g = sns.FacetGrid(df, row='method', col='n_group')
g.map(sns.scatterplot, 'd', 'p_val', s=10)
g.set_titles(template='{row_name}, {col_name} groups')
for ax in g.axes.flatten():
ax.tick_params(length=5, labelleft=True, labelbottom=True)
ax.axhline(0.05, c=plt.cm.tab10(3), lw=1)
g.fig.tight_layout(w_pad=1, h_pad=2)
g.fig.show()
g.fig.savefig('tmp.png')
追記
各グループのデータ点が少ない場合
各グループのデータ点を100ずつから10ずつに変えた場合も調べました。違いが分かりやすいようdの最大値も1から2に変えました。Tukey法とBonferroni法の差が少しだけ広がりましたが、大きな差はないままでした。
各グループの平均が非等間隔の場合
各グループの平均が非等間隔の場合も調べました。上のコードで
Y = X + np.linspace(0, d, n_group)
の箇所を
Y = X + d * np.linspace(0, 1, n_group) ** 0.5
としました。結論は変わりませんでした。
最小のp値の比較
最小のp値を比較したところ、BH法の値 (=Bonferroni法の値) よりわずかにTukey法の値の方が小さいことが確認できました。前述のとおり、2番目以降はBH法の方が拾いやすいです。
-
WelchのT検定やU検定の後にBoferroni法を使うという意味です。 ↩