【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

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

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）

著者情報： 斎藤 勇哉

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