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


1. 目的
2. 準備
2.1. data.csv
3. ソースコード
4. 結果
4.1. 対応なし検定の結果（independent_result.csv）


1. 目的

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

2. 準備

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

pip3 install scipy
pip3 install statsmodels

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

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

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

3. ソースコード

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使った方が、きれいになったやろうな。。。

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
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()

4. 結果

コードを実行すると、対応あり検定の場合、「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）を以下に示す。

4.1. 対応なし検定の結果（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

著者情報： 斎藤 勇哉

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

Print Friendly, PDF & Email
Python, Statistics, 未分類
2022-09-05リンクURL コメントを残す

コメントを残す

This site uses Akismet to reduce spam. Learn how your comment data is processed.