【初心者向け】基礎&実践プログラミング

初心者がつまづきやすいところ、最短で実力が身につく方法をお伝えします。

【Python】対応のあり・なし、パラメトリック・ノンパラメトリック検定と多重比較補正

f:id:AIProgrammer:20210115193208p:plain

動かしながら学ぶ PyTorchプログラミング入門

動かしながら学ぶ PyTorchプログラミング入門

  • 作者:斎藤勇哉
  • 発売日: 2020/11/30
  • メディア: Kindle版

目次

目的

  • Pythonを使って、対応のあり・なし、パラメトリック・ノンパラメトリック検定と多重比較補正

準備

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

pip3 install scipy
pip3 install statsmodels

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

labelは、2群分だけつけることができ、0, 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

ソースコード

data.csvが用意出来たら、検定の種類(test_type) を、対応あり(paired)か対応なし(independent )かを選択する。

検定の手法は、以下を採用している。

  • 正規性の検定として、 Kolmogorov-Smirnov testおよびShapiro-Wilk test
  • パラメトリック(para)検定では、Student's t-test
  • ノンパラメトリック(nonpara)検定として、Wilcoxon signed-rank test(対応あり)とMann–Whitney U test(対応なし)
  • 多重比較補正として、Bonferroniと Benjamini/Hochberg

コードは、Dictionary使った方が、きれいになったやろうな。。。

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


def statistical_test(data, test_type):
    # Define
    group0 = data[data['label'] == 0].drop('label', axis=1)
    group1 = data[data['label'] == 1].drop('label', axis=1)

    # Kolmogorov-Smirnov test
    kstest_group0, kstest_group1 = stats.kstest(group0, 'norm'), stats.kstest(group1, 'norm')
    # Shapiro-Wilk test
    shapiro_group0, shapiro_group1 = stats.shapiro(group0), stats.shapiro(group1)

    if test_type == 'independent':
        # Independent test
        para = stats.ttest_ind(group0, group1)  # Independent Student's t-test
        nonpara = stats.mannwhitneyu(group0, group1, alternative='two-sided')  # Mann–Whitney U test
        para_pval = para.pvalue[0]  # extract value from array
    else:
        # Paired test
        para = stats.ttest_rel(group0, group1)  # Paired Student's t-test
        nonpara = stats.wilcoxon(group0, group1)  # Wilcoxon signed-rank test
        para_pval = para.pvalue
    return kstest_group0.pvalue, kstest_group1.pvalue, shapiro_group0.pvalue, shapiro_group1.pvalue, para_pval, nonpara.pvalue

def correct_pval(pval):
    bonferroni = statsmodels.stats.multitest.multipletests(pval, alpha=0.05, method='bonferroni', is_sorted=False, returnsorted=False)  # Bonferroni
    fdr_bh = statsmodels.stats.multitest.multipletests(pval, alpha=0.05, method='fdr_bh', is_sorted=False, returnsorted=False)  # Benjamini/Hochberg (non-negative)
    return bonferroni[1], fdr_bh[1]

def main():
    # Define
    test_type = 'independent'  # independent or paired
    kstest_group0 = []
    kstest_group1 = []
    shapiro_group0 = []
    shapiro_group1 = []
    para_pval = []
    nonpara_pval = []
    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']]
        temp_kstest_group0_pval, temp_kstest_group1_pval, temp_shapiro_group0_pval, temp_shapiro_group1_pval, temp_para_pval, temp_nonpara_pval = statistical_test(data, test_type)
        kstest_group0.append(temp_kstest_group0_pval)
        kstest_group1.append(temp_kstest_group1_pval)
        shapiro_group0.append(temp_shapiro_group0_pval)
        shapiro_group1.append(temp_shapiro_group1_pval)
        para_pval.append(temp_para_pval)
        nonpara_pval.append(temp_nonpara_pval)

    # Correct pval
    para_bonferoni_pval, para_fdr_bh_pval = correct_pval(para_pval)
    nonpara_bonferoni_pval, nonpara_fdr_bh_pval = correct_pval(nonpara_pval)
        
    # Save result
    result_df = pd.DataFrame(
        {"Variable": input_df.drop(['ID', 'label'], axis=1).columns, 
        'KS group0': kstest_group0, 
        'KS group1': kstest_group1, 
        'SW group0': shapiro_group0, 
        'SW group1': shapiro_group1, 
        'para uncorrected':para_pval, 
        'para Bonferroni':para_bonferoni_pval, 
        'para FDR BH':para_fdr_bh_pval, 
        'nonpara uncorrected':nonpara_pval, 
        'nonpara Bonferroni':nonpara_bonferoni_pval, 
        'nonpara FDR BH':nonpara_fdr_bh_pval})
    result_df.to_csv("{}_result.csv".format(test_type), index=False)


if __name__ == "__main__":
    main()

結果

コードを実行すると、対応あり検定の場合、「paired_result.csv」が出力され、対応なし検定の場合、「independent_result.csv」が出力される。

各列の情報は、次の通り。

  • Variable : 検定した変数
  • KS group0 : label=0にけるKolmogorov-Smirnov testのP値
  • KS group1 : label=1にけるKolmogorov-Smirnov testのP値
  • SW group0 : label=0にけるShapiro-Wilk testのP値
  • SW group1 : label=1にけるShapiro-Wilk testのP値
  • para uncorrected :パラメトリック検定のP値
  • para Bonferroni : Bonferroni補正したP値 (FWE補正)
  • para FDR BH : Benjamini/Hochberg補正したP値 (FDR補正)
  • nonpara uncorrected : ノンパラメトリック検定のP値
  • nonpara Bonferroni : Bonferroni補正したP値 (FWE補正)
  • nonpara FDR BH : Benjamini/Hochberg補正したP値 (FDR補正)

例として、対応なし検定の結果(independent_result.csv)を以下に示す。

対応なし検定の結果(independent_result.csv)

Variable KS group0 KS group1 SW group0 SW group1 para uncorrected para Bonferroni para FDR BH nonpara uncorrected nonpara Bonferroni nonpara FDR BH
variable1 6.55E-10 5.15E-16 0.070608 0.76393 0.44383 1 0.44383 0.504789 1 0.504789
variable2 8.10E-10 8.03E-16 0.627943 0.136468 0.018309 0.091544 0.045772 0.017217 0.086084 0.043042
variable3 5.16E-12 6.71E-19 0.964636 0.867577 0.279383 1 0.349229 0.271301 1 0.339126
variable4 5.42E-12 5.06E-19 0.855986 0.480649 0.068376 0.341879 0.11396 0.078682 0.393412 0.131137
variable5 2.86E-10 1.19E-16 0.897948 0.228705 0.003365 0.016823 0.016823 0.000966 0.004831 0.004831


頑張れ!喝!!の代わりにB!ブックマークを押していただけるとただただうれしいです(^^)! ↓

動かしながら学ぶ PyTorchプログラミング入門

動かしながら学ぶ PyTorchプログラミング入門

  • 作者:斎藤勇哉
  • 発売日: 2020/11/30
  • メディア: Kindle版