多重比較におけるTukey法とBonferroni法とBH法の比較

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 |

## 結果

![tmp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2787764/38c102af-b100-e69c-487f-5a4fb176e7fc.png)

上から順に、補正なし、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]。

[^1]: WelchのT検定やU検定の後にBoferroni法を使うという意味です。

もちろん、この結果からBH法が良いと考えることも可能です。ただし、FWER制御とFDR制御の違いに注意して下さい。つまり、陽性と判定されたもののうち5%程度は誤検出かもしれないがそれでも構わない、という場合にBH法が使えます。

## ソースコード

```Python
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法の差が少しだけ広がりましたが、大きな差はないままでした。

![tmp_n10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2787764/997b435c-b3ea-6d7d-fce6-7c285fb81571.png)

### 各グループの平均が非等間隔の場合

各グループの平均が非等間隔の場合も調べました。上のコードで
```Y = X + np.linspace(0, d, n_group)```
の箇所を
```Y = X + d * np.linspace(0, 1, n_group) ** 0.5```
としました。結論は変わりませんでした。

![tmp_nonregular.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2787764/a7d629bb-28a5-1ff5-ee78-b6a8b535a19b.png)

### 最小のp値の比較

最小のp値を比較したところ、BH法の値 (=Bonferroni法の値) よりわずかにTukey法の値の方が小さいことが確認できました。前述のとおり、2番目以降はBH法の方が拾いやすいです。