[python][pandas][Statistics]pythonでの多重検定について-False Discovery rateを制御するBenjamini and Hochbergの方法を中心に解説-

Table of Contents

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

計算方法とq値

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など方法があることはあるらしい。良いやり方がわかったら記載する。