【Python】1要因で分類される対応のない多群の検定・事後検定（post hoc test）

1. 目的
2. 準備
2.1. data.csv
3. ソースコード
4. 結果
4.1. パラメトリック検定の結果（para＿result.csv）
4.2. パラメトリック検定の結果（nonpara＿result.csv）

1. 目的

• Pythonで、「1要因で分類される対応のない多群の検定（ANOVA or Kruskal-Wallis test）」と「事後検定（Tukey-Kramer or Steel-Dwass test ）」

2. 準備

統計解析には、scipyとscikit-posthocsを用いる。インストールしていない場合は、以下のコマンドでインストールする。

pip3 install scipy
pip3 install scikit-posthocs

以下のように、データを準備する。この時、必ず「ID」と「label」という列名を入れる。その他の比較する変数（variable）の列名および列数は、任意。

labelは、3群分だけつけることができ、0, 1, 2でラベル付けすること。

2.1. data.csv

3. ソースコード

data.csvが用意出来たら、検定の種類（test_type） を、パラメトリック（para）かノンパラメトリック（nonpara）かを選択する。

• パラメトリック（para）を選択した場合、pre-hocでANOVA、post-hocでTukey-Kramer testを実行。
• ノンパラメトリック（nonpara）を選択した場合、pre-hocでKruskal-Wallis test、post-hocでSteel-Dwass testを実行。

設定ができら、実行する。

いろいろいじっていたら、コードが汚くなってしまった。。。

import pandas as pd
from scipy import stats
import scikit_posthocs as sp


def statistical_test(data, test_type):
    if test_type == 'para':
        # parametric
        prehoc = stats.f_oneway(
            data[data['label'] == 0].drop('label', axis=1),
            data[data['label'] == 1].drop('label', axis=1),
            data[data['label'] == 2].drop('label', axis=1))  # ANOVA test
        posthoc = sp.posthoc_tukey(
            data,
            val_col=data.drop('label', axis=1).columns[0],
            group_col='label')  # Turkey test
    else:
        # nonpara metric
        prehoc = stats.kruskal(
            data[data['label'] == 0].drop('label', axis=1),
            data[data['label'] == 1].drop('label', axis=1),
            data[data['label'] == 2].drop('label', axis=1))  # Kruskal-Wallis test
        posthoc = sp.posthoc_dscf(
            data,
            val_col=data.drop('label', axis=1).columns[0],
            group_col='label')  # Dwass, Steel, Critchlow and Fligner all-pairs comparison test
    return prehoc, posthoc


def main():
    # Define
    test_type = 'nonpara'  # para or nonpara
    prehoc_p = []
    posthoc_p01 = []
    posthoc_p02 = []
    posthoc_p12 = []
    input_df = pd.read_csv('data.csv')

    # Test
    print("Test type: {}metric".format(test_type))
    print("Variable:")
    for col in input_df.drop(['ID', 'label'], axis=1).columns:
        print(col)
        data = input_df.loc[:, [col, 'label']]
        prehoc, posthoc = statistical_test(data, test_type)
        prehoc_p.append(prehoc.pvalue)
        posthoc_p01.append(posthoc[0][1])
        posthoc_p02.append(posthoc[0][2])
        posthoc_p12.append(posthoc[1][2])

    # Save result
    result_df = pd.DataFrame(
        {"Variable": input_df.drop(['ID', 'label'], axis=1).columns, "0 vs 1 vs 2": prehoc_p, "0 vs 1": posthoc_p01, "0 vs 2": posthoc_p02, "1 vs 2": posthoc_p12})
    result_df.to_csv("{}_result.csv".format(test_type), index=False)


if __name__ == "__main__":
    main()

4. 結果

コードを実行すると、パラメトリック（para）検定の場合、「para＿result.csv」が出力され、ノンパラメトリック（nonpara）検定の場合、「nonpara＿result.csv」が出力される。

ラベルlabel（0, 1, 2）ごとの、差の検定結果としてp値が、比較する変数ごとにまとめられている。

4.1. パラメトリック検定の結果（para＿result.csv）

4.2. パラメトリック検定の結果（nonpara＿result.csv）

著者情報： 斎藤 勇哉

順天堂大学医学部 大学院医学研究科 放射線診断学講座所属
脳MRI 画像解析が専門であり、テーマは①神経変性疾患の機序解明、②医用人工知能の開発、③多施設データのハーモナイゼーション、④速読が脳に与える影響や学習効果、⑤SNS解析を用いたマーケティング戦略の改善。
医療分野に関わらず、自然言語処理・スクレイピング・データ分析・Web アプリ開発を得意とし、企業や他大学の研究を支援。
主な使用言語は、Python、Shell Script、MATLAB、HTML、CSS