pythonでの多重検定の方法について
調べても意外にreferenceが見つからない。最初はStataでやろうとしたが、DataFrameの操作がやりにくく、python(pandas)で簡単に解決したので記録として。
多重検定の必要性
例えば薬剤の解析などで、1,000種類の薬剤の解析を行うとする。ある2群（interventionとcontrol）で薬剤の効果に差が無いかを見る時などのシチュエーションが想定できる。この時、有意差があるかないかという基準に\(p = 0.05\)を用いると、一つ一つの検定ではtype I error（帰無仮説が実際は真であるのに棄却してしまう過誤）が起こる確率は0.05であるが、1,000回の検定を行うと、それぞれの検定が独立している場合、False positiveが出現する確率は\(1-(1-0.05)^{10000} \sim 1\)となり、ほとんど確実にFalse positiveが出現してしまう。これを補正するために、様々な方法がある。もっともよく知られているものにはBonferroni補正があり、これは有意水準を（上の例で言うと）\(0.05/1000\)にしてしまうというものである。これはこれで良いが、conservativeというか、True positiveが除かれてしまうという問題がある。
False Discovery Rate (FDR)
そこでFalse Discovery Rate（FDR）を制御して、False positiveが少々混ざってもいいので、True positiveを積極的に見つけようというのがBenjamini and Hochbergの方法である。具体的にやってみると非常に簡単である。まずはDataFrameを作る。
import pandas as pd
dt = pd.DataFrame({"ID":[0, 1, 2, 3, 4, 5, 6, 7],
"Outcome1":[0.6, 0.4, 0.11, 0.99, 0.34, 0.41, 0.05, 0.87],
"Outcome2":[0.32, 0.44, 0.65, 0.34, 0.21, 0.45, 0.55, 0.5],
"Outcome3":[0.01, 0.99, 0.88, 0, 0.89, 0.67, 0.78, 0.3],
"Outcome4":[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7],
"Outcome5":[0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2],
"Category":[1, 0, 0, 1, 0, 0, 0, 1]
})
dt.set_index("ID")これは次のようなDataFrameとなる。
単純にそれぞれのOutcomeに対してt検定を行うと、次のようになる。
from scipy import stats
_, p1 = stats.ttest_ind(dt[dt["Category"] == 0]["Outcome1"], dt[dt["Category"] == 1]["Outcome1"], equal_var = False)
_, p2 = stats.ttest_ind(dt[dt["Category"] == 0]["Outcome2"], dt[dt["Category"] == 1]["Outcome2"], equal_var = False)
_, p3 = stats.ttest_ind(dt[dt["Category"] == 0]["Outcome3"], dt[dt["Category"] == 1]["Outcome3"], equal_var = False)
_, p4 = stats.ttest_ind(dt[dt["Category"] == 0]["Outcome4"], dt[dt["Category"] == 1]["Outcome4"], equal_var = False)
_, p5 = stats.ttest_ind(dt[dt["Category"] == 0]["Outcome5"], dt[dt["Category"] == 1]["Outcome5"], equal_var = False)それぞれの結果は下。
#p1 0.01767066192081987
#p2 0.45954205828448846
#p3 0.005538249034800617
#p4 0.9126960857932551
#p5 0.9126960857932551\(p < 0.05\)を有意水準にとると、Outcome1とOutcome3は有意差あり、ということになる。Benjamini and Hochbergの方法では、まずp値を昇順に並べる。
#p3 0.005538249034800617
#p1 0.01767066192081987
#p2 0.45954205828448846
#p4 0.9126960857932551
#p5 0.9126960857932551そして、施行の回数をnとして、上から順に\(n/i\)を掛けていく（本当は下から考えるのだが、便宜的にこう説明する）。
| p値 | 計算方法 | q値 |
| p3 | \(p3 \times 5/1\) | 0.08835330960409934 |
| p1 | \(p1 \times 5/2\) | 1.148855145711221 |
| p2 | \(p2 \times 5/3\) | 0.009230415058001028 |
| p4 | \(p4 \times 5/4\) | 1.140870107241569 |
| p5 | \(p5 \times 5/5\) | 0.9126960857932552 |
False positiveをどの程度許容するかの基準をq値と呼んで、p値と区別する。q値の有意水準を0.05とすると、有意と考えていたOutcome3は実は有意ではないという判断になる。q値の有意水準としては、0.1を用いることも多いので、結局結論として変わらない、という判断になることもある。pythonのcodeとしては、下のようになる。
pvalue_list = []
#ttestし、結果をpvalue_listに格納
for i, j in enumerate(dt.columns[0:-1]):
_, temp = stats.ttest_ind(dt[dt["Category"] == 0][j], dt[dt["Category"] == 1][j], equal_var = False)
if np.isnan(temp):
continue
else:
pvalue_list.append([temp, j])
#小さい順に並び替え
pvalue_list = sorted(pvalue_list, key = lambda x:x[0], reverse = False)
#Benjamini and Hochbergのq valueを添加
for i, j in enumerate(pvalue_list):
j.append(j[0]*(len(pvalue_list)/(i+1)))
#sliceのためDataFrameに変換
pvalue_list = pd.DataFrame(pvalue_list)
pvalue_list.rename(columns = {0:"p", 1:"Name", 2:"q"}, inplace = True)
pvalue_list
#q < 0.1を満たすものを抽出
pvalue_list[pvalue_list["q"] < 0.1]
#q < 0.2を満たすものを抽出
pvalue_list[pvalue_list["q"] < 0.2]
#q < 0.3を満たすものを抽出
pvalue_list[pvalue_list["q"] < 0.3]dt.columnsとしてDataFrameを横方向に眺めていくことで、sequentialに解析を行い、順にp値を格納していくイメージ。多重検定に対する補正法の問題については、色々あるがここでは触れない。
以下のリンクを参考にした。
https://sci-tech.ksc.kwansei.ac.jp/~tohhiro/dataScience/dataScience3.pdf
Stataではあまりうまくできなかったが、multprocなど方法があることはあるらしい。良いやり方がわかったら記載する。
コメント