【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

ID	label	variable1	variable2	variable3	variable4	variable5
HC01	0	0.414136	0.414828	0.59054	0.595829	0.46384
HC02	0	0.420073	0.418989	0.557235	0.569612	0.403623
HC03	0	0.4426	0.427471	0.553696	0.534986	0.399849
HC04	0	0.425205	0.437847	0.556361	0.561195	0.443079
HC05	0	0.436437	0.414882	0.573477	0.567248	0.442456
HC06	0	0.422437	0.412239	0.569644	0.574168	0.40977
HC07	0	0.446981	0.43914	0.581699	0.600997	0.451804
HC08	0	0.373986	0.398213	0.551188	0.552817	0.397256
HC09	0	0.446453	0.43022	0.588915	0.595184	0.415407
HC10	0	0.442333	0.41936	0.600606	0.582284	0.381209
HC11	0	0.42958	0.422916	0.567845	0.566097	0.412061
HC12	0	0.415686	0.418727	0.618557	0.616877	0.477297
HC13	0	0.390931	0.39321	0.559449	0.567954	0.414736
HC14	0	0.384482	0.389655	0.501869	0.509973	0.382196
HC15	0	0.39486	0.39924	0.513154	0.53733	0.346855
HC16	0	0.440336	0.4355	0.58567	0.587052	0.396427
HC17	0	0.42659	0.425409	0.582757	0.573467	0.42602
HC18	0	0.377684	0.409419	0.544228	0.539868	0.398466
HC19	0	0.400032	0.411624	0.557889	0.571171	0.429983
HC20	0	0.42478	0.396172	0.530866	0.52576	0.446319
HC21	0	0.380697	0.410013	0.533047	0.528432	0.41215
MS01	1	0.443497	0.432588	0.607397	0.604467	0.492007
MS02	1	0.43027	0.401163	0.576346	0.582719	0.390165
MS03	1	0.444296	0.403215	0.598099	0.579954	0.375619
MS04	1	0.395553	0.402554	0.553537	0.576749	0.36633
MS05	1	0.398515	0.409528	0.560015	0.561606	0.361341
MS06	1	0.409198	0.404887	0.583462	0.560014	0.39752
MS07	1	0.462465	0.452396	0.562167	0.546603	0.495371
MS08	1	0.427022	0.417257	0.559981	0.554595	0.414634
MS11	1	0.430838	0.439219	0.575243	0.574212	0.393967
MS12	1	0.393287	0.353302	0.546746	0.554518	0.345266
MS13	1	0.377094	0.383602	0.533615	0.521044	0.308438
MS14	1	0.340822	0.326417	0.522482	0.508222	0.265355
MS15	1	0.387435	0.396106	0.559567	0.568845	0.374467
MS16	1	0.409211	0.398715	0.524179	0.522349	0.406039
MS18	1	0.454714	0.424418	0.550706	0.537401	0.396089
MS19	1	0.419654	0.409647	0.550946	0.57421	0.435959
MS20	1	0.430882	0.430659	0.543696	0.550128	0.408815
MS23	1	0.410008	0.406569	0.537766	0.541704	0.357089
MS24	1	0.397079	0.386071	0.581993	0.533163	0.372379
MS25	1	0.396569	0.384555	0.542021	0.542999	0.380674
MS28	1	0.426207	0.406263	0.565803	0.551525	0.356858
MS29	1	0.415905	0.403991	0.564604	0.526406	0.357118
MS30	1	0.409348	0.396887	0.531576	0.530988	0.353894
MS39	1	0.361565	0.36709	0.521476	0.535813	0.27698
MS40	1	0.450216	0.427933	0.606333	0.613349	0.468908
MS41	1	0.396214	0.375978	0.529719	0.527544	0.331372
MS42	1	0.435776	0.420972	0.574889	0.572943	0.418034
MS43	1	0.42056	0.414211	0.529785	0.526468	0.401055
MS46	1	0.393243	0.398004	0.554618	0.564712	0.363159
MS48	1	0.407397	0.411551	0.589607	0.571527	0.418004
MS50	1	0.37233	0.383161	0.504406	0.519802	0.282472
MS51	1	0.357431	0.336885	0.477739	0.511818	0.23742
MS53	1	0.433408	0.413104	0.575848	0.556432	0.386467
MS54	1	0.411719	0.377535	0.545174	0.53249	0.356906
NMO01	2	0.432394	0.428111	0.583013	0.583175	0.430922
NMO02	2	0.398525	0.41198	0.500598	0.504638	0.385263
NMO03	2	0.365173	0.357349	0.551819	0.56269	0.3884
NMO04	2	0.42515	0.407857	0.572765	0.554362	0.430229
NMO05	2	0.388122	0.401515	0.557498	0.577619	0.400303
NMO06	2	0.341607	0.374971	0.507936	0.525418	0.39888
NMO07	2	0.42467	0.391724	0.581764	0.57302	0.381552
NMO08	2	0.441589	0.415995	0.549388	0.583331	0.413541
NMO09	2	0.434015	0.450014	0.564149	0.572999	0.41399
NMO10	2	0.395299	0.403833	0.572748	0.58591	0.407776
NMO11	2	0.360151	0.366118	0.492654	0.498297	0.341672
NMO12	2	0.423888	0.422382	0.575465	0.577919	0.443287
NMO13	2	0.287746	0.315613	0.478305	0.504997	0.236532
NMO14	2	0.404791	0.414589	0.578303	0.567639	0.409896
NMO15	2	0.43772	0.412525	0.562516	0.569668	0.43291
NMO16	2	0.459309	0.444111	0.588605	0.58474	0.437424
NMO17	2	0.400355	0.402395	0.544468	0.545933	0.373918
NMO18	2	0.377378	0.391985	0.522774	0.5222	0.350065

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）

Variable	0 vs 1 vs 2	0 vs 1	0 vs 2	1 vs 2
variable1	0.27281861	0.766279	0.248509	0.491958
variable2	0.07479037	0.07835	0.177403	0.9
variable3	0.35186936	0.552365	0.33802	0.816177
variable4	0.20801274	0.181545	0.533445	0.846405
variable5	0.01154528	0.008206	0.297811	0.412675

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

Variable	0 vs 1 vs 2	0 vs 1	0 vs 2	1 vs 2
variable1	0.529886	0.757013	0.49883	0.81342
variable2	0.055597	0.044349	0.243408	0.9
variable3	0.467262	0.50926	0.547232	0.9
variable4	0.194727	0.180888	0.756737	0.593363
variable5	0.003333	0.002696	0.337092	0.20822

著者情報：斎藤勇哉

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