「入門はじめての統計解析」をPythonで写経 Vol.17 ～ ５章「はじめてのノンパラメトリック検定」④アンサリー・ブラッドレイ検定、２×２クロス集計表、コーエンのκ統計量

５章「はじめてのノンパラメトリック検定」

書籍の著者　石村貞夫 先生

この記事は、書籍「入門はじめての統計解析」５章「はじめてのノンパラメトリック検定」の Python写経活動 を取り扱います。

書籍の図・表・計算を淡々とPython化する写経シリーズです。
５章の統計的検定は「ノンパラメトリック」。
母集団の確率分布を仮定しない統計的検定です。

この記事は５章の残りのテーマに取り組みます。

• アンサリー・ブラッドレイ検定（２標本のばらつきの検定）

• $${\boldsymbol{2 \times 2}}$$ クロス集計表と独立性の検定、イエーツの補正、フィッシャーの直接法

• コーエンの $${\kappa}$$ 統計量（評価の一致度をはかる）

ChatGPTの活用も継続してまいります！

では書籍を開いて統計解析の旅に出発です🚀

はじめに

このブログシリーズは、書籍「入門はじめての統計解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「統計解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

５章 はじめてのノンパラメトリック検定

この記事は５章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

５章で用いるライブラリをインポートします。

### インポート

# 数値計算
import collections      # python標準ライブラリ
import itertools        # python標準ライブラリ
import math             # python標準ライブラリ
import numpy as np
import pandas as pd

# 統計
import scipy.stats as stats
import pingouin as pg
from statsmodels.stats.descriptivestats import sign_test  # 符号検定

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

イントロダクション

ChatGPTが語りかける「カジュアルめのイントロ」です。

「この２つのグループ、なんかばらつき方が違う気がする…」
「カテゴリで分けたら傾向に差があるかも？」
「２人の評価、パッと見そろってるけど…ほんとに信頼していいの？」

そんな“なんとなく気になる”感覚に答えてくれるのが、今回の検定たちです。

• アンサリー・ブラッドレイ検定は、データの「ばらつき」が違うかどうかを確かめたいときに。

• $${2 \times 2}$$ クロス集計表を使った検定たちは、カテゴリどうしの「関係」や「偏り」があるかどうかを見つけてくれます。

• コーエンの $${\kappa}$$ 統計量は、「一致してるように見える評価」が、どれだけ信頼できるかを測ってくれる便利な道具です。

今回は、

• データの「ばらつき」に注目するアンサリー・ブラッドレイ検定

• 「関係」や「偏り」を見抜く $${2 \times 2}$$ クロス集計表の検定

• 「一致」の確かさを評価するコーエンの $${\kappa}$$ 統計量

この３つの視点で、データともうちょっと深く仲良くなってみましょう😊

Section 5.7 その他のノンパラメトリック検定

アンサリー・ブラッドレイ検定に取り組みます。
２つのデータの「バラツキ」に差があるかの検定です。

■ ウォーミングアップ
ChatGPTがサマリーを作ってくれました。

📝 アンサリー・ブラッドレイ検定とは？

データの平均は同じっぽい。でも、ばらつき方が違うかも？
そんなときに使えるのが、アンサリー・ブラッドレイ検定です。

🔍 どんな検定？

• 2つの独立したグループに対して

• 中央値（位置）は同じと仮定した上で

• ばらつき（スケール）に差があるかを調べる検定です。

🎯 こんなシーンで

• 実験の結果が「平均は同じなのに、安定性に違いがありそう」

• Ａグループはデータがまとまっていて、Ｂグループはバラつきが大きい

• 分散の差を比べたいけれど、正規性の仮定は避けたい

✅ 特徴ポイント

• ノンパラメトリック → 分布の形を仮定しない

• 平均ではなく → ばらつきに注目

• ウィルコクスン検定が“位置の差”を見るのに対し、
  アンサリー・ブラッドレイ検定は“広がりの差”を見ます。

✨ まとめ一言

「位置が同じなら、今度は広がりを比べてみよう」
それが、アンサリー・ブラッドレイ検定の発想です。

■ アンサリー・ブラッドレイ検定：検定統計量の公式 p.244
アンサリー・ブラッドレイ検定の検定統計量 $${A}$$ は「中心からの位置順位の和」です。
ウィルコクスンの順位和検定に似ていますがちょっと違います。

テキストの公式に代えて「順位和の仕組み」を図表を交えて書きます。
※テキストの奇数の場合と異なる可能性があることをご承知おきください。

２つのグループが次のようなデータであるとします。

２つのグループ全体で順位付けしましょう。
データを小さい順に並べて順位を付与します。

順位付けが特徴的です。
中央値に近いほど順位が大きく、中央値から遠くなるほど順位が小さくなります。

順位付けは簡単。両端から順位を設定します。
最も値の小さいものと大きいものを順位 $${1}$$ とします。
２番目に小さいものと大きいものを順位 $${2}$$ とします。

２つのグループに分かれている元の表に順位列を付け足しましょう。
そして１つのグループの順位の和をとります。

グループ１の順位の和は $${8}$$ です。
順位和が検定統計量 $${A}$$ ですので、検定統計量 $${A=8}$$ です。

◆

仮説は次のとおりです（両側検定の場合）。

なお棄却域は「アンサリー・ブラッドレイの数表」より求めます（テキスト p.294 の数表 7 参照）。
検定の都度、「数表」をめくって棄却域を取得するのは大変だと感じ…
自作関数で算出するようにしました。

■ 自作関数の留意事項
計算ロジックの決定的な情報を見つけられず、自作関数で適切な計算ができているか自信がありません…
特に次のケースは自信がないです。
・小標本の場合
・同順位（タイ）がある場合

■ アンサリー・ブラッドレイ検定統計量算出関数の実装

アンサリー・ブラッドレイ検定の検定統計量 $${A}$$ を算出する関数を定義します。
順位計算に scipy.stats を利用します。

### アンサリ・ブラッドレイ検定 p.244

# アンサリ・ブラッドレイ検定の検定統計量Aの算出関数
def stat_ansari_bradley(list1, list2):
    # 標本サイズの算出
    N1, N2 = len(list1), len(list2)
    # ２つの標本をあわせて順位付け（昇順で1からN1+N2まで） ※scipy.stats利用
    rank = stats.rankdata([*list1, *list2], method='average')
    # 最小値・最大値が1位になるように順位付けを補正
    if (N1 + N2) % 2 == 0:  # N1+N2が偶数の場合
        scores = [x if x <= (N1 + N2)/2 else (N1 + N2 + 1) - x for x in rank]
    else:                   # N1+N2が奇数の場合
        scores = [x if x <= (N1 + N2 + 1)/2 else (N1 + N2 + 1) - x for x in rank]
    # 検定統計量Aの算出: グループ1(標本1)の順位和
    A = sum(scores[:N1])
    # 戻り値: 検定統計量A, 全データの順位和、グループ1の順位和, グループ2の順位和
    return {'A': A, 'rank': rank, 'All_score': scores, 
            'Group_1_score': scores[:N1], 'Group_2_score': scores[N1:]}

続いて標本サイズが大きい場合の検定関数を定義します。
検定統計量 $${TA}$$ が標準正規分布近似することを利用しています。
放送大学さんのWebサイトの「検定統計量 $${A}$$ から検定統計量 $${TA}$$ への変換式」を参考にいたしました。
ありがとうございます！
こちらのリンクからご覧いただけます。

# 大標本の検定統計量TAの算出関数
# 標準正規分布近似でp値・棄却限界値を算出できる
# 参考サイト: 
# https://www.ouj.ac.jp/mijika/tokei/contents/sub_contents/c01_11_00.xml#link2_05
def stats_ansari_bradley_large_samples(list1, list2):
    # 検定統計量Aの算出
    result = stat_ansari_bradley(list1, list2)
    A = result['A']
    # 検定統計量TAの算出
    N1, N2 = len(list1), len(list2)
    a = N1*(N1 + N2 + 2) / 4
    b = ((N1*N2*(N1 + N2 - 2)*(N1 + N2 + 2)) / (48*(N1 + N2 - 1)))**(1/2)
    TA = (A - a) / b
    # 戻り値: 検定統計量A, 検定統計量TA
    return {'A_stat': A, 'TA_stat': TA}

# アンサリ・ブラッドレイ検定（大標本用）の実行関数
def ansari_bradley_test_large_samples(list1, list2, alpha=0.05,
                                      alternative='two-sided'):
    # 設定と準備
    std_norm_dist = stats.norm(loc=0, scale=1)  # 標準正規分布の設定

    # 検定統計量A, TA, TAの標準正規分布の上側確率の算出 ※scipy.stats利用
    A, TA = stats_ansari_bradley_large_samples(list1, list2).values()
    prob = std_norm_dist.sf(x=abs(TA))
    
    # 両側検定・片側検定の別に、棄却限界値c_val, p値p_valの算出
    match alternative:
        case 'two-sided':
            c_val = std_norm_dist.ppf(q=1 - alpha/2)
            c_val = -c_val, c_val
            p_val = prob * 2
        case 'greater':
            c_val = std_norm_dist.ppf(q=1 - alpha)
            p_val = prob
        case 'less':
            c_val = std_norm_dist.ppf(q=alpha)
            p_val = prob
    
    # 戻り値: 検定統計量A, TA, 棄却限界値, 有意水準, p値, 両側・片側検定
    return {'A_stat': A, 'TA_stat': TA, 'c_value': c_val,
            'alpha': alpha, 'p_value': p_val, 'alternative': alternative}

■ アンサリー・ブラッドレイ検定の実行【標本サイズが大きい場合】
まず仮想データを作成します。

### 標本サイズが大きい場合のアンサリー・ブラッドレイ検定 ※標本サイズ 130

# 仮想データの作成

# 設定: 乱数生成器
rng = np.random.default_rng(seed=5) # 有意性なし 1, 3, 6, 7, 8

# データの作成
list1 = rng.gamma(shape=8, scale=0.7, size=70) + 1.6
list2 = rng.gamma(shape=9, scale=0.8, size=60)

# 位置（中央値）の確認
print(f'G1の中央値 {np.median(list1):.2f}, G2の中央値 {np.median(list2):.2f}')

【実行結果】
中央値が近い２グループのデータです。

データを可視化して、直感的に「バラツキの違い」を見ましょう。

# ヒストグラムで可視化
sns.histplot(data=dict(グループ1=list1, グループ2=list2), bins=22, ec='white')
plt.xlabel('時間（秒）', fontsize=12);

【実行結果】
グループ１（青）のほうがバラツキが小さい感じがいたします。

両側検定でアンサリー・ブラッドレイ検定を実行します。

# アンサリ・ブラッドレイ検定 大標本の場合（正規分布近似）
result = ansari_bradley_test_large_samples(list1, list2, alternative='two-sided')
result

【実行結果】
有意水準 $${\alpha=5\%}$$ で有意はと言えず、帰無仮説を棄却できません。
バラツキが異なるとは言えないです。

出力内容を簡単に説明します。

• A_stat：検定統計量 $${A}$$

• TA_stat：検定統計量 $${TA}$$

• c_value：棄却限界値（$${TA}$$を使う）

• alpha：有意水準 $${\alpha}$$

• p_value：$${p}$$ 値
  帰無仮説のもとで検定統計量が実現値になる確率
  $${p}$$ 値が有意水準以下（または未満）のとき、帰無仮説を棄却できる

• alternative：検定方法

scipy.stats の結果と比べましょう。

# scipy.stats利用 ※タイを含む場合、ワーニング発生
# 標本サイズが共に55未満で同値がない場合，与えられたp値は正確なものとなり，
# そうでない場合はp値の正規近似が使用されます

stats.ansari(list1, list2, alternative='two-sided')

【実行結果】
検定統計量と $${p}$$ 値です。
検算結果は一致しました。

検定統計量 $${TA}$$、標準正規分布、棄却域を可視化します。

### 検定結果の可視化

## 設定と準備
ends = abs(result['TA_stat']) + 1
ends = 3 if ends < 3 else ends
left, right = -ends, ends                   # x軸の描画範囲
color = 'tab:blue'                          # 基本の色
c_val_low, c_val_high = result['c_value']   # 棄却限界値の下側・上側
TA = result['TA_stat']                      # 検定統計量TA
x_val_all = np.linspace(left, right, 1001)  # x軸の値:標準正規分布の確率密度関数用
x_val_low = np.linspace(left, c_val_low)    #        棄却域（下側）用
x_val_high = np.linspace(c_val_high, right) #        棄却域（上側）用
std_norm_dist = stats.norm(loc=0, scale=1)  # 標準正規分布の設定

## 描画
plt.figure(figsize=(7, 3))
plt.plot(x_val_all, std_norm_dist.pdf(x_val_all), label='標準正規分布')
plt.fill_between(x_val_low, 0, std_norm_dist.pdf(x_val_low), color=color,
                 alpha=0.2, label=f"有意水準{result['alpha']:.1%}棄却域")
plt.fill_between(x_val_high, 0, std_norm_dist.pdf(x_val_high), color=color,
                 alpha=0.2)
plt.axvline(TA, color='tab:red', ls='--', label=f'検定統計量{TA:.4f}')
plt.title(f"アンサリ・ブラッドレイ検定, $p$値{result['p_value']:.4f}")
plt.xticks([c_val_low, 0, c_val_high])
plt.legend();

【実行結果】

■ アンサリー・ブラッドレイ検定の実行【標本サイズが小さい場合】
ChatGPT ＆ Gemini と壁打ちして、小標本用の検定関数を作成しました。
ただし正直、自信がないです。。。
いくつかのテストデータについて scipy の ansari() の結果と比べてほぼ合っていたので、暫定公開します。
同順位（タイ）のある場合、変な感じになります。

### アンサリ・ブラッドレイ検定（小標本用）の実行関数 ★結果に自信がないです…

def ansari_bradley_test_small_samples(
        list1, list2, alpha=0.05, alternative='two-sided'):
    
    ## グループ１の検定統計量Aと全データの順位データ（中心スコア方式）の取得
    result = stat_ansari_bradley(list1, list2)
    obs_stat = result['A']        # グループ１の検定統計量A
    scores = result['All_score']  # 中心スコア方式による順位和[1, 2, …, 2, 1]
    
    ## 検定統計量Aの分布・確率の算出
    # 取りうる検定統計量Aの算出：可能な組み合わせより検定統計量Aを算出
    all_combos = list(itertools.combinations(range(len(scores)), len(list1)))
    all_stats = [sum(scores[i] for i in combo) for combo in all_combos]
    # 取りうる検定統計量Aの要素（順位和）ごとの度数・確率値を算出
    all_counter = dict(sorted(collections.Counter(all_stats).items()))
    all_nums, all_freqs = list(all_counter.keys()), list(all_counter.values())
    all_probs = dict(zip(all_nums, [a / sum(all_freqs) for a in all_freqs]))
    # グループ１の検定統計量Aの確率値を算出
    obs_prob = all_probs.get(obs_stat, 0)

    ## p値の算出
    match alternative:
        case 'two-sided':
            # 観測されたWと同じか、それより稀な（確率が小さい）ものを合計
            p_lower = sum(prob for a, prob in all_probs.items() if a <= obs_stat)
            p_upper = sum(prob for a, prob in all_probs.items() if a >= obs_stat)
            p_value = min(1.0, min(p_lower, p_upper) * 2)
        case 'less':
            # Wが大きいほどA群のばらつきが小さい → 片側（右側）
            p_value = sum(prob for a, prob in all_probs.items() if a >= obs_stat)
        case 'greater':
            # Wが小さいほどA群のばらつきが大きい → 片側（左側）
            p_value = sum(prob for a, prob in all_probs.items() if a <= obs_stat)
        case _:
            p_value = None

    return {
        'A_stat':obs_stat, 'alpha': alpha, 'p_value': p_value,
        'alternative': alternative, 'obs_prob': obs_prob, 'all_stats': all_stats,
        'all_probs': {
            'num': all_nums, 'freq': all_freqs, 'prob': list(all_probs.values())}, 
        }

続いて仮想データを作成します。

### 標本サイズが大きい場合のアンサリー・ブラッドレイ検定 ※標本サイズ 10

# 仮想データの作成
list1 = [1.1, 1.2, 1.4]
list2 = [1.0, 1.5, 1.6, 1.15, 1.25, 0.8, 1.9]

# 位置（中央値の確認）
print(f'G1の中央値 {np.median(list1)}, G2の中央値 {np.median(list2)}')

【実行結果】
中央値が近い２グループのデータです。

データを可視化して、直感的に「バラツキの違い」を見ましょう。

# ヒストグラムで可視化
sns.histplot(data=dict(グループ1=list1, グループ2=list2), bins=22, ec='white')
plt.xlabel('時間（秒）', fontsize=12);

【実行結果】
グループ１（青）のほうがバラツキが小さい感じがいたします。

両側検定でアンサリー・ブラッドレイ検定を実行します。

# アンサリー・ブラッドレイ検定の実行

# 検定の設定：両側検定
alternative = 'two-sided'

# 検定の実行
result = ansari_bradley_test_small_samples(list1, list2, alternative=alternative)

# 検定統計量とp値の表示
print('statistic =', result['A_stat'], ', p_value =', result['p_value'])

【実行結果】
検定統計量 $${A}$$ と $${p}$$ 値を表示しています。
有意水準 $${\alpha=5\%}$$ で有意はと言えず、帰無仮説を棄却できません。
バラツキが異なるとは言えないです。

scipy.stats の結果と比べましょう。

# scipy.stats 利用
stats.ansari(list1, list2, alternative=alternative)

【実行結果】
検定統計量と $${p}$$ 値です。
検算結果は一致しました。

標本サイズが大きい場合の関数を使ってみましょう。

# 大標本用の検定を試す
ansari_bradley_test_large_samples(list1, list2, alternative=alternative)

【実行結果】
$${p}$$ 値の値は異なるものの、有意水準 $${5\%}$$ で有意にならないことは一緒でした。

検定統計量 $${TA}$$の分布、棄却域を可視化します。
ChatGPT に原案をつくってもらいました。
こちらも自信がないです。

# 検定統計量 A の帰無分布グラフ描画 ★自信がないです…
# 引数 検定結果 res
def plot_ansari_exact_distribution_with_correct_rejection(res):

    ## 設定と準備
    nums = res['all_probs']['num']    # 取りうる検定統計量Aの値（順位和）
    probs = res['all_probs']['prob']  # 取りうる検定統計量Aの確率
    a_pairs = list(zip(nums, probs))  # 上記の値と確率のペアリスト
    obs_stat = res['A_stat']          # グループ１の検定統計量A
    alt = res['alternative']          # 検定方法（両側検定など）
    alpha = res['alpha']              # 有意水準α
    rejection = []                    # 棄却域となる値（順位和）を格納するリスト

    ## 棄却域の判定と検定方法（文字列）の収録
    # 両側検定の場合
    if alt == "two-sided":
        ## 棄却域 rejection の算出
        # 確率の小さい順にペアリストを並び替え
        a_pairs_sorted = sorted(a_pairs, key=lambda x: x[1])
        cum_prob = 0
        # 有意水準αに達するまで確率の小さい値（順位和）を取得
        for w, p in a_pairs_sorted:
            if cum_prob + p <= alpha:
                rejection.append(w)
                cum_prob += p
            else:
                break
        # タイトルの文字列の設定
        method_name = '両側検定'
    # 片側検定（下側）の場合
    elif alt == "less":
        ## 棄却域 rejection の算出
        # 確率の小さい順にペアリストを並び替え
        a_pairs_sorted = sorted(a_pairs, key=lambda x: x[0])
        cum_prob = 0
        # 有意水準αに達するまで確率の小さい値（順位和）を取得
        for w, p in a_pairs_sorted:
            if cum_prob + p <= alpha:
                rejection.append(w)
                cum_prob += p
            else:
                break
        # タイトルの文字列の設定
        method_name = '片側検定（下側）'
    # 片側検定（上側）の場合
    elif alt == "greater":
        ## 棄却域 rejection の算出
        # 確率の大きい順にペアリストを並び替え
        a_pairs_sorted = sorted(a_pairs, key=lambda x: x[0], reverse=True)
        cum_prob = 0
        # 有意水準αに達するまで確率の小さい値（順位和）を取得
        for w, p in a_pairs_sorted:
            if cum_prob + p <= alpha:
                rejection.append(w)
                cum_prob += p
            else:
                break
        # タイトルの文字列の設定
        method_name = '片側検定（上側）'

    # 棄却域用のリストから重複する値を取り除く
    rejection_set = set(rejection)

    ## 描画処理
    # 描画領域の設定
    fig, ax = plt.subplots(figsize=(8, 5))
    # バーの色の設定（棄却域：赤、通常：青）
    colors = ['tab:red' if n in rejection_set else 'tab:blue' for n in nums]
    # 確率の棒グラフの描画
    bars = ax.bar(nums, probs, color=colors, alpha=0.7)
    # 棒グラフの上に確率値（％）を表示
    ax.bar_label(bars, fmt=lambda x : '{:.{}f}%'.format(x*100, 2))
    # グループ１の検定統計量Aを示す垂直点線の描画
    ax.axvline(obs_stat, color='black', ls='--', lw=2,
               label=f'検定統計量 $A$ = {obs_stat}')
    # 修飾
    ax.set_xlabel("検定統計量 A", fontsize=12)
    ax.set_ylabel("確率", fontsize=12)
    ax.set_title('アンサリー・ブラッドレイ検定の帰無分布\n' 
                 f'{method_name}  $p$ 値 {res["p_value"]:.4f}')
    ax.legend()
    plt.tight_layout()

    return fig, ax

## アンサリー・ブラッドレイ検定の帰無分布の可視化
plot_ansari_exact_distribution_with_correct_rejection(res);

【実行結果】
赤いバーが棄却域です。
検定統計量を示す垂直点線が棄却域に含まれていない様子が分かります。

ハラハラ・ドキドキしながら関数を作りました。
とても楽しく制作時間を過ごせたので、ＯＫといたしましょう🍀

２×２クロス集計表と独立性の検定

テキスト p.237（スピアマンの順位相関係数による検定）にひっそり佇むテーマです！
セル数が少ない $${2 \times 2}$$ クロス集計表（２行２列）の独立性検定に関する「あるあるテーマ」です。

■ ウォーミングアップ
ChatGPTが「独立性の検定」「イエーツの補正」「フィッシャーの直接法」のサマリーを作ってくれました。

🧪 「独立性の検定」「イエーツの補正」「フィッシャーの正確確率検定」とは？

✅ 1. カイ二乗検定（独立性の検定）

✅ 2. イエーツの補正つきカイ二乗検定

✅ 3. フィッシャーの正確確率検定（Fisher’s Exact Test）

🧩 3つの検定のちがいまとめ

🧭 どう使い分けたらいいの？

✨ まとめ

それぞれの特性を知っておけば、データの量や性質に応じて、安心して検定が選べるようになります。
“使い分ける力”が、統計の頼もしさを一段と引き出してくれますよ😊

■ ３つの検定の実践

データ $${X, Y}$$ の $${2 \times 2}$$ クロス集計表のイメージから始めましょう。
次の例は、２つの部署の社員 $${18}$$ 人がアンケート「AIを活用しているか」に回答した「はい」「いいえ」を集計したものです。

次の図は、集計セルを「変数（記号）」に置き換えて抽象化したものです。

ではデータの作成から始めましょう。

### 仮想データの作成：生成AI活用アンケートの集計結果
data1 = np.array([[6, 1],
                  [3, 8]])
pd.DataFrame(data1, index=['第１営業部', '第２営業部'],
             columns=['活用している', '活用していない'])

【実行結果】
２つの部署でAI活用の有無が違って見えます。
部署とAI活用の有無は関連がある（独立でない）かもしれません。

🔢 独立性の検定
テキストの検定統計量の数式をお借りして、上表の記号で書き換えます。

独立性の検定の検定統計量

は、自由度 $${1}$$ の $${\chi^2}$$ 分布に従います。

scipy.stats の chi2_contingency で独立性の検定を実行しましょう。
帰無仮説は「２つの部署のAI活用有無は独立である（関連がない）」、対立仮説は「２つの部署のAI活用有無は関連がある」です。

# 独立性の検定 scipy.stats利用

# 戻り値をprintする関数の定義
def print_res(res):
    print('statistic:', res.statistic)
    print('pvalue   :', res.pvalue)
    print('dof      :', res.dof)
    print('expected_freq:')
    print(res.expected_freq)

res = stats.chi2_contingency(data1, correction=False)
print_res(res)

【実行結果】
$${p}$$ 値は $${<0.05}$$ なので有意水準 $${5\%}$$ で有意となり、帰無仮説は棄却されます。
２つの部署とAI活用有無は関連がありそうです。
営業第１部の方が活用している感じですね！

🔢 イエーツの補正を行うときの独立性の検定
テキストの検定統計量の数式をお借りして、上表の記号で書き換えます。

イエーツの補正を行うときの独立性の検定の検定統計量は

となります。

分母の $${\tfrac{N}{2}}$$ は各セルの計算値から $${0.5}$$ を差し引く操作です。
この「差し引き」を用いて検定統計量の値が小さくすることで、標本サイズが小さい場合に生じる可能性がある「検定統計量が過大になること」、「結果として $${p}$$ 値が過小になること」を抑えられるようです。

scipy.stats で「イエーツの補正」を行う独立性の検定を実行しましょう。
引数 correction=True と設定します。

# イエーツの補正
res = stats.chi2_contingency(data1, correction=True)
print_res(res)

【実行結果】
なんと！
$${p}$$ 値は $${>0.05}$$ なので有意水準 $${5\%}$$ で有意とは言えず、帰無仮説を棄却できません！
イエーツの補正を行わない場合と異なる結論になりました。。。

確かにイエーツの補正を行わないと検定統計量が過大・$${p}$$ 値が過小になるのかもしれません。

そこで！

🔢 フィッシャーの正確確率検定
テキストが「フィッシャーの直接法」と呼ぶ検定です。
統計WEBさんの情報によると、次の場合に該当する場合、独立性の検定結果が正確でないことがあるそうです。

仮想データを点検しますと…
「４つのセルのうち２つのセル、つまり 50% のセルの期待度数が５未満」になっていました！
※独立性の検定の出力結果の「expected_freq」で期待度数を確認できます。

テキストの数式をお借りして、上表の記号で書き換えます。
この数式は セル a （１行１列目のセル）の値に関する超幾何分布の確率値を求めるもののようです。

$${!}$$ は階乗の記号です。
例えば $${3!}$$ は $${3 \times 2 \times 1 = 6}$$ です。

scipy.stats の fisher_exact でフィッシャーの正確確率検定を実行しましょう。

# フィッシャーの正確確率検定
res = stats.fisher_exact(data1, alternative='two-sided')
print('statistic:', res.statistic)
print('pvalue   :', res.pvalue)

【実行結果】
なんと！
$${p}$$ 値は $${<0.05}$$ なので有意水準 $${5\%}$$ で有意となり、帰無仮説は棄却されます！

いったいどの検定を信じればいいのか…
（検定前に検定手法を決めておきましょう）

テキストの確率値の公式と $${p}$$ 値の関係を深堀りしましょう。
セル a の確率値を求めます。

# テキストの公式でセルaの確率値を計算

# 階乗関数の定義
fc = lambda x: math.factorial(x)

# フィッシャーの直接法関数の定義
fisher_exact_formula = lambda a, b, c, d: (
    fc(a+b) * fc(c+d) * fc(a+c) * fc(b+d)
    / (fc(a+b+c+d) * fc(a) * fc(b) * fc(c) * fc(d))
)

# セルaの確率値を算出
fisher_exact_formula(*data1.flatten())

【実行結果】
セル a = 6 の確率値は $${0.02376}$$ です。

仮想データの場合に「セル a の取りうる値」に関する超幾何分布の確率を計算します。

# 仮想データの超幾何分布の確率を計算

## 設定と準備
# 超幾何分布のパラメータの設定
M = data1.sum()
n = data1[0].sum()
N = data1[:, 0].sum()
# 仮想データを用いた超幾何分布の取りうる値の算出（下端・上端）
start, end = stats.hypergeom.support(M, n, N)

## 超幾何分布の確率質量関数の算出 
hg_pmfs = stats.hypergeom.pmf(np.arange(start, end+1), M, n, N)
pd.DataFrame(dict(値=range(start, end+1), 確率=hg_pmfs))

【実行結果】
セル a の取りうる値と確率の関係が分かりました。
a = 6 の確率値は、テキストの公式の計算結果と同じ $${0.02376}$$ になりました。

いよいよ、a = 6 の確率値と $${p}$$ 値の関係に迫ります！
両側検定のケースで調べます。
実は、上の「取りうる」表のうち、a = 6 の確率値 $${0.02376}$$ 以下の確率値を合算すると $${p}$$ 値になるのです。
次のコードで確かめましょう。

# scipy.statsのp値計算を追う（両側検定）

# scipy.stats のp値を取得
res = stats.fisher_exact(data1, alternative='two-sided')
p_value = res.pvalue
print(f'scipyのfisher_exactのp値: {p_value:.8f}')
print('-'*50)

# 超幾何分布の確率質量関数より、セル「a」の確率以下になる確率値を列挙
under_cell_a = hg_pmfs[hg_pmfs <= hg_pmfs[data1[0, 0]]]
print(f'セルaの確率値 {hg_pmfs[data1[0, 0]]:.8f} 以下の確率値:\n{under_cell_a}')
print(f'これらの確率値の合計 = p値: {sum(under_cell_a):.8f}')

【実行結果】
scipy の $${p}$$ 値とテキストの確率計算式が繋がりましたね！

フィッシャーの正確確率検定は $${2 \times 2}$$ クロス集計表（のセル a の値）が従う超幾何分布に基づいて「正確に」確率計算するものなのです。

コーエンの κ 統計量

テキスト p.243（ケンドールの順位相関係数による検定）の片隅でずっと待っていたテーマです！
コーエンの $${\kappa}$$ 統計量は、２人が評価者が同一サンプルを評価するときの「評価の一致度」をはかる統計量です。

■ ウォーミングアップ
「初めてコーエンの $${\kappa}$$（カッパ）統計量にふれるあなた向け」へと、ChatGPTがプレゼントしてくれた

🎯 ざっくり！コーエンのκ統計量のここがすごい！

ぜひお読みください！

✅ 1. 「一致率」を、ただの割合で終わらせない！

たとえば、２人の評価者が「合ってたかどうか」を見るとき、
「80％一致した！」って言われたら、なんとなく良さそうに見えますよね。

でもちょっと待って。
偶然でも一致することってあるんじゃない？🤔

たとえば、「はい／いいえ」が半々の質問を、２人ともなんとなく「はい」に偏って答えたら、偶然一致が増えるかもしれない。

🔽

✅ 2. “偶然の一致”を引き算して、純粋な一致度を測る！

ここで登場するのが コーエンの $${\boldsymbol{\kappa}}$$ 統計量（kappa）！
なんとこの指標…

「実際の一致率」から「偶然一致するはずの割合」を引いて、
“本当の意味で、ちゃんと一致した分だけ”を評価してくれるんです✨

🔽

✅ 3. だから「信頼できる合意かどうか」がわかる！

たとえば：

• 医師Ａと医師Ｂの診断がどれくらい一致しているか

• 教師２人の成績評価の合意度はどれくらいか

• ２つの分類モデルの判定はどれくらい揃っているか

こういった場面で、「一致率が高い＝いい」とは限らない問題を、
$${\kappa}$$ 統計量はスマートに解決してくれます🧠💡

🔽

📊 どれくらいの値なら“よい一致”なの？

🌟 まとめ：コーエンの $${\boldsymbol{\kappa}}$$ のここがすごい！

• ✅ ただの一致率よりもずっと信頼できる！

• ✅ 偶然の一致を除いた「純粋な合意度」を測れる！

• ✅ 評価・診断・分類の信頼性を定量化できる！

■ コーエンの $${\kappa}$$ 統計量の実践
コーエンの $${\kappa}$$ 統計量は、クロス集計表の「一致度」に関する指標です。
※統計的検定ではありません。

データのクロス集計表のイメージから始めましょう。
次の例は、新種フルーツ $${25}$$ 個の評価サンプルについて、２人の評価者がグレード $${A, B, C}$$ で評価した結果をクロス集計したものです。

次の図は「変数（記号）」に置き換えて抽象化したものです。
テキストの表5.6.3「コーエンの $${\kappa}$$ 統計量」にほぼ合わせています。
※標本サイズを $${N}$$ に変更しました。

テキストのコーエンの $${\kappa}$$ 統計量の公式をお借りします。

よく見かける公式は、観測一致率 $${P_0}$$、偶然一致率 $${P_e}$$ を用いて 次のように表されます。

テキストの公式との関連は

であり、一般的な公式の分母・分子に $${N^2}$$ を掛けるとテキストの公式に変形できます。

■ 例題でコーエンの $${\kappa}$$ 統計量に近づく

では追加インポートから始めてまいりましょう！

# 追加インポート
from sklearn.metrics import cohen_kappa_score
from statsmodels.stats.inter_rater import cohens_kappa

続いて評価データを設定して、クロス集計表（データフレーム）を作成します。

## 仮想データの作成

# X氏、Y氏の評価データの作成
data_x = ['A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 
          'B', 'B', 'B', 'B', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C']
data_y = ['A', 'A', 'A', 'B', 'A', 'C', 'C', 'B', 'A', 'A', 'B', 'B', 'B',
          'C', 'C', 'C', 'C', 'A', 'A', 'B', 'C', 'C', 'C', 'C', 'C']

# クロス集計表の作成
data2 = pd.crosstab(data_y, data_x, rownames=['Y'], colnames=['X'])
data2

【実行結果】
数値を見ても一致度は分かりません。。。

一致率を計算しましょう。

# 一致率の算出
np.mean([x == y for x, y in zip(data_x, data_y)])

【実行結果】
$${52\%}$$ が一致しています。
半分程度の一致率は一致していると言えるのでしょうか…

ちなみに評価が「一致」するのは、左上からの対角線上のセルです。
$${25}$$ 個の評価サンプルのうち、$${13}$$ 個の評価が一致しています。

🔢 自作関数で $${\kappa}$$ 統計量を計算
最初にテキストの公式で計算します。
関数の定義と計算実行を行います。

## コーエンのκ統計量の算出 テキストの公式利用

# コーエンのκ統計量算出関数の定義
# 引数 クロス集計表のnumpy配列
def calc_cohen_kappa(cross_tab):
    # 計算要素の算出
    N = cross_tab.sum(axis=None)  # 全データの個数
    c = cross_tab.sum(axis=0)     # 縦計：列単位のデータの個数
    r = cross_tab.sum(axis=1)     # 横計：行単位のデータの個数
    # コーエンのκ統計量の算出
    kappa = (N * sum(np.diag(cross_tab)) - sum(c * r)) / (N**2 - sum(c * r))
    
    # 戻り値：コーエンのκ統計量
    return kappa

# 仮想データをnumpy配列に変換
data2_np = data2.values

# コーエンのκ統計量の算出
calc_cohen_kappa(data2_np)

【実行結果】
$${\kappa = 0.2874}$$ です。
「Landis & Koch」の分類による目安では「やや一致している」です。
弱めの一致なのでしょうか。

Python のライブラリでコーエンの $${\kappa}$$ 統計量を計算しましょう。

🔢 Python ライブラリで $${\kappa}$$ 統計量を計算
ここでは、scikit-learn と statsmodels を利用します。

🖲️scikit-learn の cohen_kappa_score()
引数は２つの評価データです。

## コーエンのκ統計量の算出 scikit-learn利用
kappa = cohen_kappa_score(data_x, data_y)
print(kappa)

【実行結果】
自作関数の結果と一致しました！

🖲️statsmodels の cohens_kappa()
引数はクロス集計表（データフレーム）です。

## コーエンのκ統計量の算出 statsmodels利用
result = cohens_kappa(data2)
print(result.kappa)

【実行結果】
自作関数の結果と一致しました！

statsmodels の結果 result にはさまざまな値が格納されています。

result

【実行結果】

記事の最後をChatGPTに締めくくってもらいましょう！
詩人版です。

📘 ChatGPTのひとこと：

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。