見出し画像

「入門はじめての統計解析」をPythonで写経 Vol.16 ~ 5章「はじめてのノンパラメトリック検定」③スピアマン・ケンドールの順位相関係数による検定

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

書籍の著者 石村貞夫 先生


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

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

この記事は スピアマンの順位相関係数による検定 と ケンドールの順位相関係数による検定 に取り組みます。
ChatGPTの活用も継続してまいります!

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

ノンオイルドレッシングのイラスト:「いらすとや」さんより

はじめに


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

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

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


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

5.5 スピアマンの順位相関係数による検定
5.6 ケンドールの順位相関係数による検定

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

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

### インポート

# 数値計算
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つのグループに差があるかどうか」をやさしく確かめるために、平均に頼らず、順位や符号の視点から比べる検定を学びました。

今回のテーマは、「2つの変数のあいだに、関係があるかどうか」。
つまり、“差”ではなく、“つながり”を見つける検定です。

「値が増えるにつれて、もう一方も増えている気がする」
「この2つ、なんとなく一緒に動いてるように見える」

そんな直感を、順位の並びから確かめるのが、
スピアマンの順位相関係数
ケンドールの順位相関係数
です。

値の大小そのものではなく、順位のパターンに注目して関係性を測るこの2つの検定は、ノンパラメトリックならではのやわらかさで、データの“つながり”を教えてくれます。

今回も、自作関数とPythonライブラリを活用しながら、“比べる”から“一緒に変わる”へ──そんな視点の広がりを感じてもらえたら嬉しいです。


「スピアマンの順位相関係数」と「ケンドールの順位相関係数」は記事 Vol.2 の「Section 1.4 いろいろな順位相関係数」で実践しています。
相関係数の計算を対象にしております(検定ではありません)。
よかったら一緒にお読みください!

Section 5.5 スピアマンの順位相関係数による検定

■ ウォーミングアップ
ChatGPTがスピアマンの順位相関係数による検定のサマリーを作ってくれました。


🧭 スピアマンの順位相関係数による検定とは?

2つの変数の関係を調べたいとき、ピアソンの相関係数はよく使われますが、「データが正規分布じゃないかも」「直線的な関係じゃなさそう」というときは注意が必要です。

そんなときに使えるのが、スピアマンの順位相関係数です。

この検定は、元のデータを順位に変換してから、相関を調べる方法です。
値の大きさそのものではなく、“並び方”の傾向を見ることで、数値のかたちにとらわれずに「一緒に増えてるかな?減ってるかな?」を確かめられます。

実際の分析では、アンケートのスコアや順序尺度、非線形な関係など、平均や分布に頼らず「並びのパターンから関係を探したい」ときにぴったりです。


■ スピアマンの順位相関係数による検定の仕組み
スピアマンの順位相関係数による検定はざっくり、パラメトリック検定「ピアソンの積率相関係数による無相関の検定」と対応するノンパラメトリック検定です。

スピアマンの順位相関係数による検定は次の2点の特徴を持っています。

  • 対応のある2変数について、変数ごとに順位データに変換し、「スピアマンの順位相関係数」を算出します。
    ※テキストはスピアマンの順位相関係数を検定統計量にしています。

  • スピアマンの順位相関係数を無相関の検定の検定統計量に変換して、無相関の検定を実施します。
    ※無相関の検定を用いる手順は、テキストと異なります。

■ スピアマンの順位相関係数による検定:検定統計量の公式 p.226
スピアマンの順位相関係数の公式をテキストをお借りします。

スピアマンの順位相関係数の公式は次のとおりです。

$$
r_s = 1 - \cfrac{6 \sum_{i=1}^N (a_i - b_i)^2}{N (N^2-1)}
$$

テキストの数式を引用

$${N}$$ は標本サイズ、$${a_i, b_i}$$ は対応のある2つの変数の順位データです。
テキストはスピアマンの順位相関係数 $${r_s}$$ を検定統計量とし、「スピアマンの順位相関検定の数表」を用いて棄却域を求めます(テキスト p.296 の数表 8, 9 参照)。

この記事では「無相関の検定」を用いて棄却域を求めます。
無相関の検定の検定統計量 $${T}$$ の算出には、テキスト p.190 の公式を用います。

検定統計量

$$
T(r_s, N) = \cfrac{r_s \sqrt{N-2}}{\sqrt{1 - r_s^2}}
$$

テキストの数式を一部改変して引用

は、自由度 $${N-2}$$ の $${t}$$ 分布に従います。
$${r_s}$$ にスピアマンの順位相関係数を使い、$${N}$$ は標本サイズです。

仮説は次のとおりです。

帰無仮説 $${H_0}$$:2変数 $${A, B}$$ に相関はない
対立仮説 $${H_1}$$:2変数 $${A, B}$$ に相関がある

テキストの記述を一部改変して引用

■ スピアマンの順位相関係数による検定関数の実装

公式やWebサイトの情報に基づいて「スピアマンの順位相関係数による検定関数」を定義します。
次のWebサイトにお世話になっています。
ありがとうございます!

計算の一部で scipy.stats を利用します。

### スピアマンの順位相関係数関数 p.52, 236
# 参考サイト 統計WEB: https://bellcurve.jp/statistics/course/26041.html

# スピアマンの順位相関係数の算出関数
def calc_spearman_rank_corr(list1, list2):

    ## 設定と準備
    # 標本サイズ: 2つの標本サイズは同じ
    N = len(list1)

    ## 順位データの作成 scipy.stats利用 ※順位付けは昇順・降順のどちらでもOK
    rank1 = stats.rankdata(list1)
    rank2 = stats.rankdata(list2)

    ## 順位相関係数の算出
    # タイ(同順位)の件数カウント
    tx = [v for v in collections.Counter(rank1).values() if v > 1]
    ty = [v for v in collections.Counter(rank2).values() if v > 1]
    
    # タイが無い場合の順位相関係数の算出
    if (tx == []) & (ty == []):
        # テキストの数式に従って、右辺の分母, rsを算出
        numerator = 6 * sum((x - y)**2 for x, y in zip(rank1, rank2))
        # 戻り値: スピアマンの順位相関係数rs
        return 1 - numerator / (N * (N**2 - 1))
    
    # タイがある場合の順位相関係数の算出: 統計WEBの数式を利用
    else:
        # データx、データyのいずれかにタイが無い場合、0のリストを設定
        tx = tx if tx != [] else [0]
        ty = ty if ty != [] else [0]
        # 統計WEBの数式に従ってTx, Ty ,rsの右辺の分母, rs(戻り値)を算出
        Tx = (N * (N**2 - 1) - sum([t * (t**2 - 1) for t in tx])) / 12
        Ty = (N * (N**2 - 1) - sum([t * (t**2 - 1) for t in ty])) / 12
        numerator = Tx + Ty - sum((x - y)**2 for x, y in zip(rank1, rank2))
        # 戻り値: スピアマンの順位相関係数rs
        return numerator / (2 * (Tx * Ty)**(1/2))

# 無相関のt検定統計量の算出関数 p.190 無相関の検定の公式利用
def t_stat_non_corr(list1, list2):
    # 標本サイズ: 2つの標本サイズは同じ
    N = len(list1)
    # スピアマンの標本相関係数の算出
    r = calc_spearman_rank_corr(list1, list2)
    # t検定統計量の算出
    tau = r * ((N - 2) / (1 - r**2))**(1/2)
    # 戻り値: t検定統計量, 標本相関係数r
    return tau, r

# 無相関のt検定関数 ※確率計算はscipy.stats利用
def non_corr_ttest(list1, list2, alpha=0.05, alternative='two-sided'):
    # 標本サイズ ※2つの標本サイズは同じ
    N = len(list1)
    # t検定統計量と標本相関係数の算出
    t_val, r = t_stat_non_corr(list1, list2)
    # 自由度(N-2)のt分布の設定 scipy.stats利用
    t_dist = stats.t(df=N - 2)
    # 棄却限界値c_valとp値p_valの算出
    match alternative:
        case 'two-sided':
            c_val = t_dist.ppf(q=1 - alpha/2)
            c_val = -c_val, c_val
            p_val = t_dist.sf(x=abs(t_val)) * 2
        case 'less':
            c_val = t_dist.ppf(q=alpha)
            p_val = t_dist.cdf(x=t_val)
        case 'greater':
            c_val = t_dist.ppf(q=1 - alpha)
            p_val = t_dist.sf(x=t_val)
    return {'r': r, 't_value': t_val, 'c_value': c_val, 'alpha': alpha,
            'p_value': p_val}

次のサイトのデータを引用して、自作関数の確認を行います。
ありがとうございます!

# テスト
# データ引用元: https://istat.co.jp/sk_commentary/correlation-test/Spearmans-Test
      
# データ
list1 = [3, 3, 3, 3, 4, 2 ,4, 4, 2, 5]
list2 = [4, 3, 2, 2, 2 ,3, 4, 4, 4, 5]

# スピアマンの順位相関係数による検定の実行
non_corr_ttest(list1, list2, alternative='greater')

【実行結果】
$${p}$$ 値が $${>0.05}$$ なので、有意水準 $${5\%}$$で有意とは言えず、帰無仮説を棄却できません。

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

  • r:スピアマンの順位相関係数=検定統計量 $${r_s}$$

  • t_value:無相関の検定で使う $${t}$$ 値

  • c_value:無相関の検定の棄却限界値

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

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

scipy.stats で検算します。

# scipy.stats で答え合わせ ※無相関の検定
stats.spearmanr(list1, list2, alternative='greater')

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

■ スピアマンの順位相関係数による検定(両側検定)例題 p.238
p.238 例)の2つの教育方法をはかる試験の成績データが相関するかどうか、スピアマンの順位相関係数による検定で調べます。

仮説は次のとおりです(両側検定)。
・帰無仮説 $${H_0}$$ 「方法1と方法2に相関はない」
・対立仮説 $${H_1}$$ 「方法1と方法2に相関がある」

テキストのデータを引用して、ヒストグラムを描画します。

### スピアマンの順位相関係数による検定 試験の成績 p.238

# データ
list1 = [5, 8, 10, 3, 4, 2, 8, 9, 7, 4]
list2 = [3, 8, 8, 4, 5, 3 ,6, 10, 8, 4]

# ヒストグラムの描画
sns.histplot(dict(演繹的方法=list1, 帰納的方法=list2), ec='white')
plt.xlabel('分数計算の試験の成績 [10点満点]', fontsize=12)
plt.ylabel('人数', fontsize=12);

【実行結果】
重ねたヒストグラムでは高得点の部分で相違があるものの、その他の部分はほぼ同じ値になっています。
関連があるように見えます!

① 自作関数で検定実行
有意水準 $${5\%}$$、片側検定(下側)でスピアマンの順位相関係数による検定を実行します。

# 関数利用

# スピアマンの順位相関係数による検定の実行
non_corr_ttest(list1, list2)

【実行結果】
検定統計量 $${r_s=0.8230}$$ に無相関の検定を適用すると、棄却限界値$${2.3060}$$ より大きい(テキストによると棄却域に入る)ので、有意水準 $${5\%}$$ で帰無仮説は棄却されます。

つまり、方法1と方法2は 相関があると言える、です。

② Python ライブラリで検定実行
ここでは、scipy.stats と pingouin を利用します。
いずれも引数に2つのサンプルデータを渡しています。

🖲️scipy.stats

# scipy.stats利用 ※無相関の検定
stats.spearmanr(list1, list2)

【実行結果】
検定統計量は自作関数と一致しました!
$${p}$$ 値は $${<0.05}$$ であり、有意水準 $${5\%}$$ で帰無仮説は棄却されます。
方法1と方法2は 相関があると言える、です。

🖲️pingouin

# pingouin利用
pg.corr(list1, list2, method='spearman')

【実行結果】
検定統計量(スピアマンの順位相関係数)r と $${p}$$ 値は自作関数と一致しました!
$${p}$$ 値は $${<0.05}$$ であり、有意水準 $${5\%}$$ で帰無仮説は棄却されます。

Section 5.5 ケンドールの順位相関係数による検定

■ ウォーミングアップ
ChatGPTがケンドールの順位相関係数による検定のサマリーを作ってくれました。


🧭 ケンドールの順位相関係数による検定とは?

スピアマンと同じく、「順位」に注目して2つの変数の関係を調べる検定です。
ケンドールの検定では、すべてのデータペアを比べて

  • 「片方が大きければもう片方も大きい(=一致ペア)」

  • 「片方が大きいけど、もう片方は小さい(=不一致ペア)」
    のどちらになるかを1つずつ数えていきます。

その結果、一致と不一致のバランスから、2つの変数が同じ方向に動いているのかどうか(=“つながり”)を測るのです。
特にサンプル数が少ないときや、外れ値の影響を減らしたいときに、ケンドールの順位相関係数は安定して判断できるとされています。

「細かな変動よりも、全体の傾向をじっくり見たい」──
そんな場面で使いたくなる、やさしい検定です。

🧩 2つの検定のちがいは?


■ ケンドールの順位相関係数による検定の仕組み
ケンドールの順位相関係数による検定はざっくり、パラメトリック検定「ピアソンの積率相関係数による無相関の検定」と対応するノンパラメトリック検定です。
※スピアマンの順位相関係数による検定と似たポジショニング!

ケンドールの順位相関係数による検定を次のステップで検討します。
各ステップは①テキストの方法、②この記事の方法でやりかたが分かれます。

  • 対応のある2変数の2つの組について、順位の大小が同じ向き・異なる向きの個数を数え上げて、以下の方法で検定統計量を算出する

    • ①テキスト:総得点で検定統計量を算出

    • ②この記事:ケンドールの順位相関係数 $${\tau}$$ から検定統計量を算出

  • 検定統計量が従う分布・確率の算出には少なくとも2系統ある

    • ①テキスト:数表を用いて棄却限界値を取得

    • ②この記事:以下の計算式で算出

      • 標本サイズが小さい場合、順位の大小の組み合わせの度数から確率を算出する

      • 標本サイズが大きい場合、標準正規分布近似で確率を算出する

「順位の大小の向き」を簡単な例で見てみましょう。
テキストの「総得点で検定統計量を算出」をやってみます。

次の表は、対応関係のある変数 $${A, B}$$の順位データです。

4つの $${A,B}$$ ペアどうしの組み合わせは次の6通りになります。

例えば1つ目の組み合わせは $${(1, 1)}$$ と $${(2, 3)}$$ です。
変数 $${A}$$ の大小関係は $${1 < 2}$$ なので「Aの大小」は「<」です。
同様に、変数 $${B}$$ は $${1 < 3}$$、「Bの大小」は「<」です。
この組み合わせの大小の向きは両方とも「<」で一致しているので、+1点を与えて、$${P}$$ が1つカウントアップされます。

4つ目の組み合わせは大小の向きが反対になるので、-1点を与えて、$${Q}$$ が1つカウントアップされます。

すべての組み合わせの手続きが終わると、$${P}$$ の合計が $${5}$$、$${Q}$$ の合計が $${1}$$ となり、検定統計量 $${S=P-Q=4}$$ が求まります。

■ ケンドールの順位相関係数による検定:検定統計量の公式 p.241
検定統計量の公式をテキストよりお借りします。

2つの組 $${(a_i, b_i), (a_j, b_j)}$$ を取り出したとき

$$
a_i < a_j,\ b_i < b_j \quad または \quad a_i > a_j,\ b_i > b_j
$$

のように、順位の大小が一致すれば、$${+1}$$ 点を与える。
逆に

$$
a_i < a_j,\ b_i > b_j \quad または \quad a_i > a_j,\ b_i < b_j
$$

のように、順位の大小が反対になっているときには、$${-1}$$ 点を与える。
このとき

$$
P= +1 点の総数, \quad Q= -1 点の総数
$$

とおくと

$$
S = P - Q
$$

が検定統計量となります。

◆ 引用ここまで ◆

ケンドールの順位相関係数 $${\tau}$$ を用いて検定統計量 $${S}$$ を算出することもできます。

$$
S = \cfrac{N(N-1)}{2}\ \tau
$$

テキストの数式を引用

仮説は次のとおりです(両側検定)。

帰無仮説 $${H_0}$$:2変数 $${A, B}$$ は互いに独立である
対立仮説 $${H_1}$$:2変数 $${A, B}$$ は関連がある

テキストの記述を一部改変して引用

棄却域は「ケンドールの順位相関検定の数表」より求めます(テキスト p.297 の数表 11 参照)。

以上がテキストの検定統計量 $${S}$$ および数表を検定に用いる場合の概要です。

■ ケンドールの順位相関係数による検定関数の実装

この記事は以下の方法でケンドールの順位相関係数による検定関数を実装します。
テキストの方法と異なっていますのでご注意下さい。

  • ケンドールの順位相関係数 $${\tau}$$ から検定統計量を算出

  • 検定統計量が従う分布・確率を以下の計算式で算出

    • 標本サイズが小さい場合、順位の大小の組み合わせの度数から確率を算出

    • 標本サイズが大きい場合、標準正規分布近似で確率を算出

「ケンドールの順位相関係数による検定関数」を定義します。
以下のWebサイトの情報にお世話になりました。
ありがとうございます!

・ケンドールの順位相関係数の算出:統計WEB

・ケンドールの順位相関係数 $${\tau}$$ の $${p}$$ 値算出関数 小標本用
 青木繁伸 先生のホームページ
 pdfファイルリンク:ケンドールのτb

「ケンドールの順位相関係数 $${\tau}$$ の算出関数」です。

### ケンドールの順位相関係数関数 p.55, 241
# 参考サイト 統計WEB: https://bellcurve.jp/statistics/course/26041.html

# ケンドールの順位相関係数の算出関数
def kendall_rank_corr(list1, list2):
    ## 設定
    # 同じ向きの組P, 逆の向きの組Q, xのタイTx, yのタイTyのカウンターの初期化
    P, Q, Tx, Ty = 0, 0, 0, 0
    # 比較した組み合わせの総数
    n = len(list1)       # 標本サイズ
    N = n * (n - 1) / 2  # 組み合わせの総数

    ## 順位データの作成 scipy.stats利用 ※順位付けは昇順・降順のどちらでもOK
    rank1 = stats.rankdata(list1)
    rank2 = stats.rankdata(list2)
    
    ## 同じ向きの組・逆の向きの組をカウント
    # 組み合わせの特定
    for i in range(n):
        for j in range(i+1, n):
            ## x,yともに組み合わせ間の差を算出
            x_comp = rank1[i] - rank1[j]
            y_comp = rank2[i] - rank2[j]
            
            ## 向きの判定
            # x,yの差が両方ともプラス or 両方ともマイナスのとき、同じ向きPをカウント
            if (x_comp > 0) & (y_comp > 0):
                P += 1
            elif (x_comp < 0) & (y_comp < 0):
                P += 1
            # x,yの差が一方がプラス・他方がマイナスのとき、逆の向きQをカウント
            elif (x_comp > 0) & (y_comp < 0):
                Q += 1
            elif (x_comp < 0) & (y_comp > 0):
                Q += 1
            # xがタイのとき、xのタイをカウント
            elif x_comp == 0:
                Tx += 1
            # xがタイのとき、xのタイをカウント
            elif y_comp == 0:
                Ty += 1
    
    ## ケンドールの順位相関係数τの算出
    # タイデータが無い場合
    if (Tx == 0) & (Ty == 0):
        tau = (P - Q) / N
    # タイデータがある場合
    else:
        tau = (P - Q) / ((N - Tx)**(1/2) * (N - Ty)**(1/2))

    ## 戻り値: ケンドールの順位相関係数τ
    return tau

「$${\tau}$$ に基づく小標本の $${p}$$ 値算出関数」です。
引用サイトのPythonコードを一部改変しています。
組み合わせの計算がとても高速です。

# ケンドールの順位相関係数τのp値算出関数 小標本用
# 引用元: 青木繁伸 先生 公開ソース
#         http://aoki2.si.gunma-u.ac.jp/Python/Kendall_tau_b.pdf

# 確率算出関数
def pkendall(q, n):
    # 標本サイズが大きい場合に警告表示して終了
    if n > 100:
        raise Exception('n must be less than 101')
    # ケンドールの検定統計量Sの個数のテーブルの作成
    len_t = int(n*(n - 1)/2 + 1)
    t = list(np.zeros(len_t, dtype=int))
    t[0] = 1
    for i in range(1, n):
        prev = t[:int(i*(i - 1)/2 + 1)].copy()
        for j in range(i):
            for k in range(len(prev)):
                t[j + 1 + k] = int(t[j + 1 + k]) + int(prev[k])
    # 確率の算出
    denom = sum(t)
    cum = sum([int(t[i]) for i in range(int(q) + 1)])
    # 戻り値: 確率
    return cum/denom

# p値算出関数
def calc_kendall_p_value(tau, n, alternative='two-sided'):
    # 組み合わせの数の位置の初期設定
    q = round((tau + 1) * n * (n - 1) / 4)
    # 両側検定・片側検定の別にp値を計算
    match alternative:
        case 'two-sided':  # 両側検定の場合
            q = min(q, int(n*(n - 1)/2) - q)
            pvalue = pkendall(q, n)
            pvalue = min(2 * pvalue, 1)
        case 'greater':    # 片側検定(上側)の場合
            pvalue = 1 - pkendall(q - 1, n)
        case 'less':       # 片側検定(下側)の場合
            pvalue = pkendall(q, n)
    # 戻り値: p値
    return pvalue

「$${\tau}$$ に基づく大標本の $${p}$$ 値算出関数」です。

# ケンドールの順位相関係数τのp値算出関数 大標本用
def calc_kendall_p_value_large(tau, N, alternative='two-sided'):
    # z値の算出
    z = tau / ((2*(2*N + 5)) / (9*N*(N - 1)))**(1/2)
    # p値の算出
    dist = stats.norm(loc=0, scale=1)
    match alternative:
        case 'two-sided':  # 両側検定の場合
            p_value = dist.cdf(-z) + dist.sf(z)
        case 'greater':    # 片側検定(上側)の場合
            p_value = dist.sf(z)
        case 'less':       # 片側検定(下側)の場合
            p_value = dist.cdf(-z)
    return p_value

「ケンドールの順位相関係数の無相関の検定関数」です。
検定を実施する際はこの関数を使います。

# ケンドールの順位相関係数の無相関の検定
def kendall_non_corr_test(list1, list2, thres=30, alternative='two-sided'):

    # 標本サイズの算出: 2つの標本サイズは同じ
    N = len(list1)
    
    # ケンドールの順位相関係数の算出
    tau = kendall_rank_corr(list1, list2)
    
    # 標本サイズの大きさ別にp値を算出
    if N <= thres:  # 小標本の場合
        p_val = calc_kendall_p_value(tau, N, alternative)
    else:           # 大標本の場合
        p_val = calc_kendall_p_value_large(tau, N)
    
    # 戻り値: ケンドールの順位相関係数tau, p値
    return {'tau': tau, 'p_value': p_val}

テストは省略いたします。

■ ケンドールの順位相関係数による検定(両側検定)例題 p.242
p.242 例)の2人のソムリエによるワインの順位付けデータが相関するかどうか、ケンドールの順位相関係数による検定で調べます。

仮説は次のとおりです(両側検定)。
・帰無仮説 $${H_0}$$ 「2人の判定順位は独立である」
・対立仮説 $${H_1}$$ 「2人の判定順位に関連がある」

テキストのデータを引用して、ヒストグラムを描画します。

### ケンドールの順位相関係数による検定 赤ワインの順位付け p.242

# データ
list1 = [5, 2, 4, 1, 3]
list2 = [4, 1, 5, 3, 2]

# 棒グラフで可視化
N = len(list1)
plt.bar([x - 0.21 for x in range(N)], list1, width=0.4, label='U氏', alpha=0.7)
plt.bar([x + 0.21 for x in range(N)], list2, width=0.4, label='O氏', alpha=0.7)
plt.xticks(ticks=range(N), labels=['A', 'B', 'C', 'D', 'E'])
plt.xlabel('赤ワインの種類', fontsize=12)
plt.ylabel('順位', fontsize=12)
plt.legend();

【実行結果】
2人の順位付けは何となく似ているように見えます。

① 自作関数で検定実行
有意水準 $${5\%}$$、両側検定でケンドールの順位相関係数による検定を実行します。

### ケンドールの順位相関係数による検定

# 関数利用
kendall_non_corr_test(list1, list2)

【実行結果】
ケンドールの順位相関係数 $${\tau}$$ と $${p}$$ 値を出力します。
$${p}$$ 値が $${>0.05}$$ であり、有意水準 $${5\%}$$ で有意とは言えず、帰無仮説を棄却できません。

つまり、2人の判定順位に 関連があるとは言えない、です。

② Python ライブラリで検定実行
ここでは、scipy.stats と pingouin を利用します。
いずれも引数に2つのサンプルデータを渡しています。

🖲️scipy.stats

# scipy.stats利用
# 帰無仮説「2つの標本に順位相関は無い」

tau, p_value = stats.kendalltau(list1, list2)
print(f'ケンドールの順位相関係数 = {tau:.6f} (p値: {p_value:.4f})')

【実行結果】
検定統計量は自作関数と一致しました!
$${p}$$ 値は $${>0.05}$$ であり、有意水準 $${5\%}$$ で帰無仮説を棄却できません。
つまり、2人の判定順位に 関連があるとは言えない、です。

🖲️pingouin

# pingouin利用
# 帰無仮説「2つの標本に順位相関は無い」

pg.corr(list1, list2, method='kendall')

【実行結果】
ケンドールの順位相関係数 r と $${p}$$ 値は自作関数と一致しました!
$${p}$$ 値は $${>0.05}$$ であり、有意水準 $${5\%}$$ で帰無仮説を棄却できません。

■ ケンドールの順位相関係数の数表を作成する
自作関数のコード化を試行錯誤しているときに作った副産物です。

関数定義です。
このWEBサイトの情報を参考にいたしました。
ありがとうございます!

### 参考: ケンドールの数表の作成関数

# ケンドールの数表の作成関数
def kendall_rank_corr_test_table(n):

    ## ケンドールの検定統計量Sの個数テーブルtの作成
    # 確率計算参考サイト: https://ushitora.net/archives/1958
    # [1, ..., n]の組み合わせ(重複なし)の作成
    a = list(itertools.permutations(list(range(1, n + 1))))
    # 各組み合わせ内のペアの作成
    c = [list(itertools.combinations(b, 2)) for b in a]
    # 各ペアについてケンドールの検定統計量Sの個数テーブルtの作成
    f = collections.Counter(
        [sum([1 if int(e[1] - e[0]) > 0 else -1 for e in d]) for d in c])
    S, t = list(f.keys()), list(f.values())

    ## 確率の計算
    # 確率の計算
    prob = [x / sum(t) for x in t]
    # 累積確率の作成
    prob_cumcum = [sum(prob[:i+1]) for i in range(len(t))]
    # 戻り値
    return {'P-Q': S, 't': t, '確率': prob, '累積確率': prob_cumcum}

# ケンドールの数表のデータフレーム作成関数
def make_kendall_number_table(n, alpha, return_df=False):
    ## ケンドールの数表のデータフレームを作成
    # ケンドールの数表のデータフレームを作成
    df = pd.DataFrame(kendall_rank_corr_test_table(n))
    # 確率と累積確率の列名を取得
    col1, col2 = df.columns[2:]

    ## データフレームのカスタマイズ by styler
    # 累積確率が有意水準に達しない場合、確率の小数点表示のみ設定
    if sum(df[col2] <= alpha) == 0:
        styler = df.style.format({col1: '{:.4f}', col2: '{:.4f}'})
    # 累積確率が有意水準に達する場合、数表のハイライトと確率の小数点表示を設定
    else:
        styler = (
            df
            .style
            .set_properties(
                **{'background-color': 'yellow'}, 
                subset=pd.IndexSlice[df[df[col2] <= alpha][col2].idxmax(), :])
            .format({col1: '{:.4f}', col2: '{:.4f}'})
        )
    
    # 戻り値:
    if return_df:  # dfを返す場合
        return styler, df
    else:          # dfを返さない場合
        return styler

ケンドールの順位相関係数検定の数表を作成します。
標本サイズ $${n=5}$$、有意水準 $${\alpha=0.05}$$(片側検定)で実行します。
$${\alpha}$$ が位置する下側累積確率の箇所をハイライトしています。

# ケンドールの数表の作成 ※n=>8くらいから処理時間が長くなります

# 設定
n = 5         # 標本サイズ
alpha = 0.05  # 有意水準(片側検定)

# ケンドールの数表の作成
styler, kendall_df = make_kendall_number_table(n, alpha, return_df=True)
styler

【実行結果】

確率を可視化しましょう。

# 確率の可視化
plt.figure(figsize=(10, 4))
sns.barplot(data=kendall_df, x='P-Q', y='確率', width=1, ec='white', alpha=0.7);

【実行結果】
左右対称の分布になっています。


記事の最後をChatGPTに締めくくってもらいましょう!

📘 ChatGPTのひとこと:

データどうしの関係を知ると、世界が少し読みやすくなります。
そして、読む力が育つと、関係も深まっていきます。
スピアマンとケンドールの検定で見つけた「やさしい相関」。
このブログとも、そんな関係を育んでいけたらうれしいです📘✨

今回の写経は以上です。


シリーズの記事

次の記事

前の記事

目次

ブログの紹介


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

1.のんびり統計

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

2.実験!たのしいベイズモデリング1&2をPyMC Ver.5で

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

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

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

4.楽しい写経 ベイズ・Python等

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

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

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

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

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

7.Python機械学習プログラミング実践記

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

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

いいなと思ったら応援しよう!

ネイピア DS 応援ありがとうございます。これからもがんばって記事を作成します!