【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

コメントを残す

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください