クラスカル・ウォリス検定 — ノンパラメトリックANOVAの理論と実装

3つ以上の群の間に差があるかどうかを調べたい場面は、実験科学や社会科学で頻繁に登場します。たとえば、3種類の肥料が植物の成長に与える影響を比較したり、4つの異なる学習法の効果を比較したりする場合です。

このような多群比較の標準的な手法は一元配置分散分析(one-way ANOVA)ですが、ANOVAには「各群のデータが正規分布に従い、分散が等しい」という仮定が必要です。データが順序尺度であったり、強く歪んでいたり、外れ値を含んでいたりする場合、この仮定は成り立ちません。

このような場面で使えるのがクラスカル・ウォリス検定(Kruskal-Wallis test)です。この検定はマン・ホイットニーU検定を3群以上に拡張したものであり、データの順位のみを使って群間差を検定します。

クラスカル・ウォリス検定を理解すると、以下のような応用が可能です。

  • 医学研究: 3段階以上の治療法の比較で、正規性を仮定できない場合の標準的な手法です
  • 教育研究: 学生の成績やアンケート評価(リッカート尺度)の多群比較に広く用いられます
  • 工業品質管理: 複数の製造ラインや条件間での品質指標の比較に適しています
  • 生態学: 複数の生息地における種の多様性指標の比較など、非正規データの多群比較に利用されます

本記事の内容

  • クラスカル・ウォリス検定の直感的理解
  • H統計量の数学的導出
  • 漸近分布と正確検定
  • 事後検定(Dunn検定)による多重比較
  • Pythonによる実装と一元配置ANOVAとの比較

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

順位ANOVAとしての直感

ANOVAを順位に適用する

クラスカル・ウォリス検定を最も直感的に理解する方法は、「データを順位に変換してからANOVAを行う」と考えることです。

まず、全群のデータを合わせて $N$ 個のデータに小さい方から順位 $1, 2, \dots, N$ を付けます。次に、各群の平均順位を計算し、それらが互いに離れているかどうかを評価します。

もし群間に差がなければ、各群の平均順位は全体の平均順位 $\frac{N+1}{2}$ の近くに集まるはずです。逆に、ある群のデータが全体的に大きければ、その群の平均順位は高くなり、他の群から離れます。

イメージとしては、$k$ 群のデータをプールして並べたとき、各群のデータが「混ざり合って」いれば差なし、「偏って」いれば差あり、と判断します。

具体的な例

3つの肥料 A, B, C を使って植物を育てた結果を考えましょう。

肥料A 肥料B 肥料C
4.5 6.2 8.1
3.8 5.8 7.5
5.1 7.0 9.2
4.2 6.5 8.8

12個のデータに順位を付けると:

3.8 4.2 4.5 5.1 5.8 6.2 6.5 7.0 7.5 8.1 8.8 9.2
順位 1 2 3 4 5 6 7 8 9 10 11 12
A A A A B B B B C C C C

各群の平均順位: A = 2.5, B = 6.5, C = 10.5。全体の平均順位 = 6.5。

Aの平均順位は低く(小さい値が多い)、Cの平均順位は高い(大きい値が多い)ことが明確です。クラスカル・ウォリス検定はこの「平均順位の散らばり」を定量化して、統計的に有意かどうかを判定します。

この直感を数式で厳密にしていきましょう。

H統計量の導出

記号の定義

$k$ 個の群があり、第 $i$ 群には $n_i$ 個の観測値があるとします。総データ数は $N = \sum_{i=1}^{k} n_i$ です。

全データを合わせて小さい順に順位を付け、第 $i$ 群の $j$ 番目のデータの順位を $R_{ij}$ とします。

第 $i$ 群の順位和と平均順位を次のように定義します。

$$ R_i = \sum_{j=1}^{n_i} R_{ij}, \quad \bar{R}_i = \frac{R_i}{n_i} $$

全体の平均順位は次のようになります。

$$ \bar{R} = \frac{N + 1}{2} $$

これは $1, 2, \dots, N$ の平均値です。

ANOVAのF統計量からの類推

一元配置ANOVAの $F$ 統計量は、群間変動と群内変動の比として定義されます。

$$ F = \frac{\text{群間変動} / (k-1)}{\text{群内変動} / (N-k)} $$

これと同じ構造を、元のデータの代わりに順位に対して適用します。

順位の群間変動(SS_between):

$$ SS_B = \sum_{i=1}^{k} n_i (\bar{R}_i – \bar{R})^2 $$

各群の平均順位が全体の平均順位からどれだけ離れているかを測ります。

順位の全変動(SS_total):

$$ SS_T = \sum_{i=1}^{k}\sum_{j=1}^{n_i} (R_{ij} – \bar{R})^2 $$

同順位がなければ、順位は $1, 2, \dots, N$ のちょうどの置換であるため、全変動は次のように計算できます。

$$ SS_T = \sum_{r=1}^{N}\left(r – \frac{N+1}{2}\right)^2 = \frac{N(N^2-1)}{12} $$

この計算は、$\sum_{r=1}^{N} r^2 = \frac{N(N+1)(2N+1)}{6}$ の公式を使い、$\bar{R} = \frac{N+1}{2}$ を代入して整理すると得られます。

H統計量の定義

クラスカル・ウォリスのH統計量は、次のように定義されます。

$$ \begin{equation} H = \frac{SS_B}{SS_T / (N-1)} = (N-1) \cdot \frac{\sum_{i=1}^{k} n_i (\bar{R}_i – \bar{R})^2}{\sum_{i=1}^{k}\sum_{j=1}^{n_i} (R_{ij} – \bar{R})^2} \end{equation} $$

同順位がなければ $SS_T = \frac{N(N^2-1)}{12}$ を代入して簡略化できます。

$$ H = (N-1) \cdot \frac{\sum_{i=1}^{k} n_i (\bar{R}_i – \bar{R})^2}{\frac{N(N^2-1)}{12}} $$

$\bar{R} = \frac{N+1}{2}$ を使ってさらに整理すると、よく教科書で見る形が得られます。

$$ \begin{equation} H = \frac{12}{N(N+1)} \sum_{i=1}^{k} \frac{R_i^2}{n_i} – 3(N+1) \end{equation} $$

この式の導出を確認しましょう。$\bar{R}_i = R_i / n_i$ なので、

$$ \sum_{i=1}^{k} n_i \bar{R}_i^2 = \sum_{i=1}^{k} \frac{R_i^2}{n_i} $$

また、$\sum_{i=1}^{k} n_i \bar{R}_i = \sum_{i=1}^{k} R_i = \frac{N(N+1)}{2}$(全順位の和)です。

$SS_B = \sum n_i \bar{R}_i^2 – N\bar{R}^2$ を計算すると、

$$ SS_B = \sum_{i=1}^{k} \frac{R_i^2}{n_i} – \frac{N(N+1)^2}{4} $$

$H = (N-1) \cdot SS_B / SS_T$ に $SS_T = \frac{N(N^2-1)}{12} = \frac{N(N-1)(N+1)}{12}$ を代入すると、

$$ H = (N-1) \cdot \frac{\sum \frac{R_i^2}{n_i} – \frac{N(N+1)^2}{4}}{\frac{N(N-1)(N+1)}{12}} $$

分子・分母から $(N-1)$ を消去し、$\frac{12}{N(N+1)}$ を掛けて整理すると、

$$ H = \frac{12}{N(N+1)} \sum_{i=1}^{k} \frac{R_i^2}{n_i} – 3(N+1) $$

が得られます。

同順位の補正

同順位が存在する場合、H統計量に補正が必要です。同順位グループが $g$ 個あり、第 $l$ グループに $t_l$ 個のデータがあるとします。

補正付きH統計量は次のようになります。

$$ \begin{equation} H_c = \frac{H}{1 – \frac{\sum_{l=1}^{g}(t_l^3 – t_l)}{N^3 – N}} \end{equation} $$

分母の補正項は、同順位がある場合に全変動 $SS_T$ が減少することを反映しています。同順位がなければ $t_l = 1$ であり、補正項はゼロになるため $H_c = H$ です。

H統計量がどのような分布に従うかを次に見ていきましょう。

帰無分布と検定手順

帰無仮説と対立仮説

$$ H_0: F_1 = F_2 = \cdots = F_k \quad \text{(全群の分布が同一)} $$

$$ H_1: \text{少なくとも1組の群の分布が異なる} $$

マン・ホイットニーU検定と同様に、クラスカル・ウォリス検定は分布全体の同一性を検定するものであり、厳密には「中央値が等しい」ことの検定ではないことに注意してください。ただし、各群の分布が同じ形状で位置のみ異なる場合(位置シフトモデル)には、中央値の等しさの検定と等価になります。

漸近分布

各群のサイズ $n_i$ が十分大きい場合、帰無仮説のもとでH統計量は漸近的にカイ二乗分布に従います。

$$ \begin{equation} H \xrightarrow{d} \chi^2(k-1) \quad \text{under } H_0 \end{equation} $$

自由度は $k – 1$(群の数マイナス1)です。

この漸近結果の直感的な理由は次の通りです。帰無仮説のもとで各群の平均順位 $\bar{R}_i$ は漸近的に正規分布に従い(中心極限定理)、H統計量は正規変量の二次形式として構成されています。正規変量の二次形式は $\chi^2$ 分布に従うため、H統計量も $\chi^2$ 分布に従います。

一般的な目安として、各群のサイズが5以上であれば $\chi^2$ 近似は実用的に十分な精度を持つとされています。

正確検定

各群のサイズが小さい場合は、正確な帰無分布を用いることが推奨されます。正確な帰無分布は、$N$ 個の順位を $k$ 群に割り当てるすべての方法 $\frac{N!}{n_1! n_2! \cdots n_k!}$ 通りについてH統計量を計算し、その分布を直接求めることで得られます。

しかし、$N$ が大きくなると組み合わせの数が爆発的に増えるため、並べ替え検定と同様にモンテカルロ法による近似が使われます。

検定手順のまとめ

  1. 全データを合わせて順位を付ける(同順位は中間順位)
  2. 各群の順位和 $R_i$ と平均順位 $\bar{R}_i$ を計算する
  3. H統計量を計算する(同順位があれば補正を行う)
  4. $\chi^2(k-1)$ 分布の上側確率からp値を求める
  5. p値が有意水準 $\alpha$ 以下なら帰無仮説を棄却する

H統計量が有意だった場合、「どの群とどの群の間に差があるか」を特定するための事後検定が必要です。

事後検定 — Dunn検定による多重比較

多重比較の必要性

クラスカル・ウォリス検定は全体的検定(omnibus test)であり、「少なくとも1組の群間に差がある」ことしか教えてくれません。具体的にどの群の組み合わせに差があるかを知るには、事後検定(post hoc test)が必要です。

ANOVAの事後検定としてテューキーのHSD検定があるように、クラスカル・ウォリス検定の事後検定としてはDunn検定(Dunn’s test)が最もよく使われます。

Dunn検定の統計量

Dunn検定は、各ペア $(i, j)$ について、平均順位の差を標準化した統計量を計算します。

群 $i$ と群 $j$ の比較のための検定統計量は次のようになります。

$$ \begin{equation} z_{ij} = \frac{\bar{R}_i – \bar{R}_j}{\sqrt{\frac{N(N+1)}{12}\left(\frac{1}{n_i} + \frac{1}{n_j}\right)}} \end{equation} $$

同順位がある場合は、分母の $\frac{N(N+1)}{12}$ を同順位補正付きの値に置き換えます。

帰無仮説のもとで $z_{ij}$ は漸近的に標準正規分布に従います。

多重比較の補正

$k$ 群の場合、$\binom{k}{2} = \frac{k(k-1)}{2}$ 個のペアの検定を同時に行うため、多重比較問題が発生します。個々のp値をそのまま有意水準と比較すると、偽陽性(帰無仮説を誤って棄却する確率)が膨らみます。

よく使われる補正方法には以下のものがあります。

ボンフェローニ補正: 各ペアのp値に比較回数 $m = \binom{k}{2}$ を掛けます。

$$ p_{\text{adjusted}} = \min(m \cdot p_{\text{raw}}, 1) $$

最も保守的な方法ですが、比較回数が多いと検出力が大きく低下します。

ホルム補正: ボンフェローニの改良版で、p値を小さい順にソートして段階的に補正します。ボンフェローニより検出力が高く、家族内過誤率(FWER)を制御します。

ベンジャミニ・ホッホベルグ(BH)補正: 偽発見率(FDR)を制御する方法で、比較回数が多い場合に検出力が高くなります。

事後検定の理論を理解したので、次にPythonで全体を実装しましょう。

Pythonによる実装

クラスカル・ウォリス検定のスクラッチ実装

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

def kruskal_wallis_test(groups):
    """
    クラスカル・ウォリス検定のスクラッチ実装

    Parameters
    ----------
    groups : list of arrays  各群のデータ

    Returns
    -------
    H : float  H統計量
    p_value : float  p値(χ²近似)
    group_info : dict  各群の統計量
    """
    k = len(groups)
    ns = [len(g) for g in groups]
    N = sum(ns)

    # 全データを結合して順位を付ける
    combined = np.concatenate(groups)
    ranks = stats.rankdata(combined, method="average")

    # 各群の順位和と平均順位
    idx = 0
    R_sums = []
    R_means = []
    for i, n_i in enumerate(ns):
        group_ranks = ranks[idx:idx + n_i]
        R_sums.append(np.sum(group_ranks))
        R_means.append(np.mean(group_ranks))
        idx += n_i

    R_bar = (N + 1) / 2  # 全体の平均順位

    # H統計量
    H = (12 / (N * (N + 1))) * sum(
        R_sums[i]**2 / ns[i] for i in range(k)
    ) - 3 * (N + 1)

    # 同順位補正
    _, counts = np.unique(ranks, return_counts=True)
    tie_sum = np.sum(counts**3 - counts)
    if tie_sum > 0:
        correction = 1 - tie_sum / (N**3 - N)
        H = H / correction

    # p値(χ²近似)
    p_value = 1 - stats.chi2.cdf(H, df=k - 1)

    group_info = {
        "n": ns, "rank_sums": R_sums, "rank_means": R_means,
        "overall_mean_rank": R_bar
    }

    return H, p_value, group_info

# --- テストデータ: 3つの教育方法のテストスコア ---
np.random.seed(42)
method_A = np.array([72, 68, 75, 70, 65, 78, 71, 69])
method_B = np.array([82, 85, 79, 88, 83, 80, 86, 81])
method_C = np.array([90, 87, 92, 85, 95, 88, 91, 89])

groups = [method_A, method_B, method_C]
group_names = ["Method A", "Method B", "Method C"]

# 自作実装
H, p_val, info = kruskal_wallis_test(groups)

# SciPyとの比較
H_scipy, p_scipy = stats.kruskal(*groups)

print("=" * 55)
print("クラスカル・ウォリス検定の結果")
print("=" * 55)
for i, name in enumerate(group_names):
    print(f"{name}: n={info['n'][i]}, 中央値={np.median(groups[i]):.1f}, "
          f"平均順位={info['rank_means'][i]:.2f}")
print(f"\n全体平均順位: {info['overall_mean_rank']:.2f}")
print(f"\n自作実装: H = {H:.4f}, p = {p_val:.6f}")
print(f"SciPy:    H = {H_scipy:.4f}, p = {p_scipy:.6f}")
print(f"\n有意水準5%で: {'棄却 → 少なくとも1組に差あり' if p_val < 0.05 else '棄却しない'}")

上の結果から、自作実装とSciPyの結果が一致していることが確認できます。3つの教育方法の間に統計的に有意な差があることが示されています。しかし、「どの方法とどの方法の間に差があるか」はまだわかりません。

Dunn検定の実装

import numpy as np
from scipy import stats

def dunn_test(groups, group_names=None, correction="bonferroni"):
    """
    Dunn検定(クラスカル・ウォリス検定の事後検定)

    Parameters
    ----------
    groups : list of arrays
    group_names : list of str
    correction : str  'bonferroni', 'holm', 'bh'
    """
    k = len(groups)
    ns = [len(g) for g in groups]
    N = sum(ns)
    m = k * (k - 1) // 2  # ペア数

    if group_names is None:
        group_names = [f"Group {i+1}" for i in range(k)]

    # 順位の計算
    combined = np.concatenate(groups)
    ranks = stats.rankdata(combined, method="average")

    idx = 0
    mean_ranks = []
    for n_i in ns:
        mean_ranks.append(np.mean(ranks[idx:idx + n_i]))
        idx += n_i

    # 同順位補正
    _, counts = np.unique(ranks, return_counts=True)
    tie_correction = np.sum(counts**3 - counts) / (12 * (N - 1))
    sigma2 = (N * (N + 1) / 12 - tie_correction)

    # 全ペアの比較
    results = []
    for i in range(k):
        for j in range(i + 1, k):
            diff = abs(mean_ranks[i] - mean_ranks[j])
            se = np.sqrt(sigma2 * (1/ns[i] + 1/ns[j]))
            z = diff / se
            p_raw = 2 * (1 - stats.norm.cdf(z))
            results.append({
                "pair": f"{group_names[i]} vs {group_names[j]}",
                "z": z, "p_raw": p_raw,
                "mean_rank_diff": mean_ranks[i] - mean_ranks[j]
            })

    # p値の補正
    p_raws = [r["p_raw"] for r in results]

    if correction == "bonferroni":
        p_adjusted = [min(p * m, 1.0) for p in p_raws]
    elif correction == "holm":
        sorted_idx = np.argsort(p_raws)
        p_adjusted = [0.0] * len(p_raws)
        for rank_pos, orig_idx in enumerate(sorted_idx):
            p_adjusted[orig_idx] = min(p_raws[orig_idx] * (m - rank_pos), 1.0)
        # 単調性の保証
        for rank_pos in range(1, len(sorted_idx)):
            orig_idx = sorted_idx[rank_pos]
            prev_idx = sorted_idx[rank_pos - 1]
            p_adjusted[orig_idx] = max(p_adjusted[orig_idx],
                                        p_adjusted[prev_idx])
    elif correction == "bh":
        sorted_idx = np.argsort(p_raws)[::-1]
        p_adjusted = [0.0] * len(p_raws)
        for rank_pos, orig_idx in enumerate(sorted_idx):
            bh_rank = m - rank_pos
            p_adjusted[orig_idx] = min(p_raws[orig_idx] * m / bh_rank, 1.0)
        # 単調性の保証
        for rank_pos in range(1, len(sorted_idx)):
            orig_idx = sorted_idx[rank_pos]
            prev_idx = sorted_idx[rank_pos - 1]
            p_adjusted[orig_idx] = min(p_adjusted[orig_idx],
                                        p_adjusted[prev_idx])

    for i, r in enumerate(results):
        r["p_adjusted"] = p_adjusted[i]

    return results

# 先ほどのデータでDunn検定
results = dunn_test(groups, group_names, correction="bonferroni")

print("\n" + "=" * 60)
print("Dunn検定(ボンフェローニ補正)")
print("=" * 60)
print(f"{'ペア':<25} {'z':>7} {'p(raw)':>10} {'p(adj)':>10} {'判定':>6}")
print("-" * 60)
for r in results:
    sig = "***" if r["p_adjusted"] < 0.001 else \
          "**" if r["p_adjusted"] < 0.01 else \
          "*" if r["p_adjusted"] < 0.05 else "n.s."
    print(f"{r['pair']:<25} {r['z']:>7.3f} {r['p_raw']:>10.6f} "
          f"{r['p_adjusted']:>10.6f} {sig:>6}")

この事後検定の結果から、どのペア間に統計的に有意な差があるかが明確になります。全体的なクラスカル・ウォリス検定で「差がある」と判定されても、すべてのペアに差があるとは限りません。Dunn検定により、具体的にどの群の組み合わせに差が存在するかを特定できます。

ANOVAとの検出力比較

正規分布と非正規分布で、クラスカル・ウォリス検定と一元配置ANOVAの検出力を比較します。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

np.random.seed(42)

k = 3  # 群数
n_per_group = 15
n_sim = 3000
alpha = 0.05
effect_sizes = np.linspace(0, 2.0, 12)

distributions = {
    "Normal": lambda shift, n: np.random.randn(n) + shift,
    "Exponential": lambda shift, n: np.random.exponential(1, n) + shift,
    "Mixed (contaminated)": lambda shift, n: np.where(
        np.random.rand(n) < 0.9,
        np.random.randn(n) + shift,
        np.random.randn(n) * 5 + shift  # 10%の外れ値
    ),
}

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

for ax, (dist_name, gen_func) in zip(axes, distributions.items()):
    power_kw = []
    power_anova = []

    for effect in effect_sizes:
        reject_kw = 0
        reject_anova = 0

        for _ in range(n_sim):
            # 群1: シフトなし, 群2: effect/2, 群3: effect
            g1 = gen_func(0, n_per_group)
            g2 = gen_func(effect / 2, n_per_group)
            g3 = gen_func(effect, n_per_group)

            # クラスカル・ウォリス
            _, p_kw = stats.kruskal(g1, g2, g3)
            if p_kw < alpha:
                reject_kw += 1

            # 一元配置ANOVA
            _, p_anova = stats.f_oneway(g1, g2, g3)
            if p_anova < alpha:
                reject_anova += 1

        power_kw.append(reject_kw / n_sim)
        power_anova.append(reject_anova / n_sim)

    ax.plot(effect_sizes, power_kw, "b-o", markersize=5, linewidth=1.5,
            label="Kruskal-Wallis")
    ax.plot(effect_sizes, power_anova, "r-s", markersize=5, linewidth=1.5,
            label="One-way ANOVA")
    ax.axhline(alpha, color="gray", linestyle="--", linewidth=0.8)
    ax.set_xlabel("Effect size", fontsize=11)
    ax.set_ylabel("Power", fontsize=11)
    ax.set_title(dist_name, fontsize=13)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_ylim(0, 1.05)

plt.tight_layout()
plt.savefig("kruskal_wallis_power.png", dpi=150, bbox_inches="tight")
plt.show()

この検出力比較から、以下の重要な結果が読み取れます。

  1. 正規分布(左図): ANOVAがクラスカル・ウォリス検定よりわずかに高い検出力を示しています。正規分布のもとではANOVAが最適な検定であるため、これは予想通りの結果です。ただし、差はごくわずかであり、クラスカル・ウォリス検定の検出力の損失は実用上無視できるレベルです

  2. 指数分布(中央図): 歪んだ分布では、クラスカル・ウォリス検定とANOVAの検出力はほぼ同等か、クラスカル・ウォリス検定がやや優位です。指数分布は右に裾が長いため、外れ値がANOVAの性能に影響を与え始めます

  3. 汚染正規分布(右図): 外れ値を含むデータでは、クラスカル・ウォリス検定がANOVAを明確に上回ります。10%の観測値が分散5倍の正規分布から生成されるこの設定では、平均に基づくANOVAは外れ値に引きずられて検出力が低下しますが、順位に基づくクラスカル・ウォリス検定はロバストに差を検出しています

結果の可視化 — 箱ひげ図と順位プロット

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

np.random.seed(42)
method_A = np.array([72, 68, 75, 70, 65, 78, 71, 69])
method_B = np.array([82, 85, 79, 88, 83, 80, 86, 81])
method_C = np.array([90, 87, 92, 85, 95, 88, 91, 89])
groups = [method_A, method_B, method_C]
group_names = ["Method A", "Method B", "Method C"]

# 全データの順位
combined = np.concatenate(groups)
ranks = stats.rankdata(combined, method="average")
N = len(combined)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

# 箱ひげ図
ax = axes[0]
bp = ax.boxplot(groups, labels=group_names, patch_artist=True,
                boxprops=dict(alpha=0.6))
colors = ["#4C72B0", "#DD8452", "#55A868"]
for patch, color in zip(bp["boxes"], colors):
    patch.set_facecolor(color)
# 個別データ点
for i, (g, color) in enumerate(zip(groups, colors)):
    jitter = np.random.normal(0, 0.04, len(g))
    ax.scatter(np.full(len(g), i + 1) + jitter, g, color=color,
               alpha=0.7, s=40, edgecolors="white", zorder=3)
ax.set_ylabel("Test score", fontsize=12)
ax.set_title("Raw data", fontsize=13)
ax.grid(True, alpha=0.3, axis="y")

# 順位プロット
ax = axes[1]
idx = 0
for i, (n_i, name, color) in enumerate(zip([len(g) for g in groups],
                                             group_names, colors)):
    group_ranks = ranks[idx:idx + n_i]
    mean_rank = np.mean(group_ranks)
    ax.scatter(np.full(n_i, i), group_ranks, color=color,
               alpha=0.7, s=60, edgecolors="white", zorder=3)
    ax.plot(i, mean_rank, "k_", markersize=20, markeredgewidth=3, zorder=4)
    ax.annotate(f"$\\bar{{R}}={mean_rank:.1f}$",
                xy=(i, mean_rank), xytext=(i + 0.25, mean_rank),
                fontsize=10, va="center")
    idx += n_i

ax.axhline((N + 1) / 2, color="gray", linestyle="--", linewidth=0.8,
           label=f"Overall mean rank = {(N+1)/2:.1f}")
ax.set_xticks(range(len(group_names)))
ax.set_xticklabels(group_names)
ax.set_ylabel("Rank", fontsize=12)
ax.set_title("Ranks (bars = group mean rank)", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis="y")

plt.tight_layout()
plt.savefig("kruskal_wallis_visualization.png", dpi=150, bbox_inches="tight")
plt.show()

この可視化から、以下のことが確認できます。

  1. 箱ひげ図(左図)は3群のスコア分布を示している。Method Aのスコアが最も低く、Method Cが最も高いことが視覚的に明確です。各群の分布の重なりが少ないことから、統計的に有意な差が検出されることが予想されます

  2. 順位プロット(右図)は各データ点の順位と群平均順位を示している。黒い横棒が各群の平均順位であり、これらが全体平均順位(灰色の破線)からどれだけ離れているかがH統計量の本質です。Method Aの平均順位は低く、Method Cの平均順位は高いことが明確に示されています

クラスカル・ウォリス検定の注意点

検定の前提条件

クラスカル・ウォリス検定は「分布フリー」と言われますが、以下の前提は必要です。

  1. 各群のデータは独立: 異なる群のデータ間に相関がないこと
  2. 各群内のデータは独立同一分布: 同じ群内のデータは同じ分布に従うこと
  3. 順序尺度以上: データに少なくとも大小関係が定義できること

分散や形状が群間で異なる場合、帰無仮説「分布が同一」が棄却されても、それが位置(中央値)の差なのかスケール(分散)の差なのか区別できません。

$k = 2$ の場合

群が2つ ($k = 2$) の場合、クラスカル・ウォリス検定はマン・ホイットニーU検定と等価になります。$H = z^2$(Uから得られる標準化統計量の2乗)の関係が成り立ちます。

効果量

検定のp値に加えて、効果の大きさを報告することが推奨されます。クラスカル・ウォリス検定に対応する効果量として、イータ二乗 $\eta^2_H$ が使われます。

$$ \eta^2_H = \frac{H – k + 1}{N – k} $$

この指標は0から1の範囲を取り、群間の差の大きさを表します。一般的な解釈の目安は、$0.01$ が小、$0.06$ が中、$0.14$ が大です。

まとめ

本記事では、クラスカル・ウォリス検定の理論をANOVAとの対比を通じて解説しました。

  • クラスカル・ウォリス検定は「順位に対するANOVA」 として直感的に理解できます。全データに順位を付け、各群の平均順位が全体の平均順位からどれだけ離れているかを評価します
  • H統計量 $= \frac{12}{N(N+1)}\sum \frac{R_i^2}{n_i} – 3(N+1)$ は、帰無仮説のもとで漸近的に $\chi^2(k-1)$ 分布に従います
  • Dunn検定は、H統計量が有意だった場合にどのペア間に差があるかを特定する事後検定です。多重比較補正(ボンフェローニ、ホルムなど)が必要です
  • ANOVAとの比較: 正規分布のもとでもクラスカル・ウォリス検定の検出力損失はわずかであり、外れ値や歪んだ分布ではANOVAよりも優れた検出力を示します

次のステップとして、以下の記事も参考にしてください。