この記事で何ができるようになるか
- データの正規性や等分散性も考慮した検定法の選択
- パラ・ノンパラメトリックの多重検定
- 一覧性の高いヒートマップによる結果図示
はじめに
自分は大学院で生物学を専攻しています。
Pythonは実験データをグラフで解析するのに非常に有用ですが、種々の検定（特に群間の有意差検定）を行うライブラリが少なく、結構困っていました。
解決策として
- Rでやる。
- Rの関数をPythonで実装する。
等が考えられます。（2に関してはRpy2を導入すれば一応できました）
しかし、どうしてもPythonでノンパラメトリック検定などを動かしたい！
そこで第３の手法 scilit_posthocs というライブラリを使ってノンパラメトリック検定を実装することにしました。
scikit_posthocsの実装にしか興味がありませんでしたら、一気に下の目次まで飛んで頂けると幸いです。
検定の基礎：検定手法の選択方法。
まず検定を始めるにあたって
以下の内容は次のJ-stageの記事を参考にしています
統計検定を理解せずに使っている人のために Ⅰ 、Ⅱ 、Ⅲ
これらの記事を私なりに解釈して記述しています。
ちょっと長くなりそうなのでまた時間があるときに更新しますが、要するに次の流れに沿います。
3群以上で対応のない(独立な)データ
↓
正規性の検定(Shapiro-Wilk検定、もしくはQQプロットでも..) → Non-paraへ
↓
等分散性の検定(Bartlett検定) → Non-paraへ
↓
一元配置分散分析(ANOVA) → Non-paraへ
↓
Tukey_HSD検定、Scheffe検定、Tukey検定(各群のnが同じ)、Dunnett検定(Control群との比較)
参考にした記事では、
- 正規性、等分散性の帰無仮説が棄却できない場合
- ANOVAで有意差が出なかった場合
はノンパラメトリック検定を推奨しています。
上のフローチャートでノンパラにいった場合
↓
等分散性の検定(Levene検定、Fligner検定)
↓
一元配置分散分析(Kruskal-Wallis検定)
↓
Steel-Dwass-Critchlow-Fligner(dscf)検定、Conover Iman検定
一応これらの検定も事後検定(posthoc)なので、一元配置分散分析で有意差がでないとだめなのでは？？ とも思いますが、上の記事によると分散分析をしなくても良い趣旨のことを書いています。
しかしながら、特にConover Imanの検定は第一種の過誤を起こしやすいのでKrusukal-Wallis検定で帰無仮説が棄却された場合のみに使うべきだと思います。ちなみにdscfとconoverを比べた場合、dscfの方がかなり検定が厳しくなります。例えばこのリンク先では、基本的にdscf検定のほうを勧めています。確かに自分の結果もconoverより,dscfの方がかなりp値は高くなりますね（検定が厳しい）。なのでdscfでも望む結果が出る場合は、dscfを適用する方が良いと思います。
#scikit-posthocsの実装
scikit_posthocsはscipyやstatsmodelsではカバーしていない、検定を多数カバーしているライブラリで非常に使いやすいです。公式HPがとても良くまとまっているので、是非チェックしてください。
公式HP
GitHubのレポジトリ
依存パッケージはNumpy,Pandas,scipy,stasmodels,matplotlib,seabornです。
scikit_posthocsはpipでインストールできます。
!pip install scikit_posthocs
主に使うことになるのは
- Tukey_HSD検定
- Tukey検定
- Scheffe検定
- ペアワイズt検定 （多分t検定を多重に繰り返してるので、行儀が悪い）
- Steel-Dwass検定
- Conover検定
です。
どの検定も(HSD以外)以下の感じで実行できます。
import scikit_posthoc as sp
import seaborn as sns
#タイタニックのデータをロード
df = sns.load_dataset("titanic")
#Steel-Dwass 検定
#val_colは値のカラム
#group_colは比較したい群のカラム
sp.posthoc_dscf(df,val_col="fare",group_col="class")
結果は次のようなデータフレームで帰ってきます。
表の中身はp値です。
検定方法の選択から図示まで、まとめて実装
以上の内容をまとめてgithubに置いときました。
rola-bio/stats_test
この中にあるstats_test.pyを作業用のディレクトリにダウンロードしてimportしてください。
そしてstats_test()を実行すると、
上のフローで示した通りにデータの正規性、等分散性を検定し分散分析まで自動で行います。
その後適した検定でデータを分析し、有意差の結果とデータの棒グラフを図示します。
デフォルトで選択されるのは、Tukey-HSD,Steel-Dwass,Conover検定のうちどれかです。
では実際にseaborn標準搭載のタイタニックの乗客データから、乗客の種類による運賃の差を、この関数を使って分析していきます。
import stats_test as st
import seaborn as sns
#タイタニックのデータをロード
df = sns.load_dataset("titanic")
df.head()
次にstats_test()でデータフレーム、検定したい値、どの要素でグループ分けをするか指定してください。
今回は乗客の種類を乗船地(embark_town)で分けてみました。
st.stats_test(df,val_col="fare",group_col="embark_town")
おっと〜〜？？これを実行するとエラーが出ちゃいました。
TypeError: '<' not supported between instances of 'float' and 'str'
どうやらfareかembark_townにエラー(nan)がありそうです。
group_colにintが混ざっていたり、Null値が混じってるとこのエラーが出ることがあります。
intエラーの時は
df["カラム名"] = df["カラム名"].astype(str)
で対処できます。
今回は下のようにdropnaでnanを取り除きました。
また、group_colで選択するカラム名に、カッコ（）が入っている場合もエラーが出るので注意してください。
st.stats_test(df.dropna(subset=["embark_town"]),val_col="fare",group_col="embark_town")
結果
**画面の見方**
画像上側の表記は検定の結果、ノンパラメトリックで不等分散だったことを示しています。
またKruskal-Wallis検定の結果で有意差がありそうだったので、Conover検定が自動で選択されました。
図の左手はデータをプロットし、右側に検定の結果をヒートマップで表しています。
どうやらすべての群間にp値<0.001以下の有意差があります。
Cherbourgで乗った人は有意におぼっちゃまですね、、、、
さらにグループを男女で分ける
おっと、私め！あろうことか男性をけなすような発言をしてしまいました。
Cherbourgで乗った人は有意におぼっちゃま
このデータは男女を区別していませんから、次は男女別に有意差検定をしてみましょう。
for sex in df["sex"].unique():
print("""
This result is from {}
""".format(sex))
df_query = df.query("sex =='{}'".format(sex))
st.stats_test(df_query.dropna(subset=["embark_town"]),
val_col="fare",group_col="embark_town")
というような結果になりました。
最初の結果と、乗船地の色分けが変わってしまったのでわかりづらいですね。。。
パッケージ内のsign_barplot()をいじっていただければ調整可能です。
何はともあれ、Cherbourgの乗客は男女ともに有意に金持ちのようです。（ぐぬぬ,,,）
しかし男性のSouthamptonとCherbourgの賃金の差はp値<0.01程度にまで上昇しました。
吉田麻也のおかげですかね？
ということでおしまいです。
ちなみにstats_test()にresult = Trueを渡すと、途中の検定の結果も表示します。
test = "テスト名"を渡すと、検定を自分で指定できます。
(もしくはstats_test.pyの関数,one_way_ANOVA()をいじっても簡単に変更できます)
そのうち個々の検定の実装も書くかもしれません。
詳しくはコードの中身を参照していただければ...ですが。