🧺

HCI研究者の心理統計ノート：クラスカル・ウォリス検定 / Kruskal-Wallis test

2023/01/09に公開

research

statistic

tech

HCI研究者の心理統計ノート

Introduction

クラスカル・ウォリス検定 (Kruskal-Wallis test) は、3つ以上の対応のないデータに差があるかどうかを検定する方法です。ネット上の多くの記事では無視されていますが等分散性が仮定出来ない場合には精度が低下するようなので (Zimmer, 2004)、等分散性が仮定出来る場合にのみ使用するのが良いようです。

HCI研究では、1要因参加者間計画 (between-subject design) の実験において、条件1,2,...,kで評価指標に統計的に有意な差があるかどうかを調べる際に使用します (e.g., 既存手法1 vs 既存手法2 vs 提案手法)。特に、質問紙の評定（リッカート尺度）を比較する場面などでよく使います。

分散分析とは違い、正規性の仮定は必要ありません。正規性が仮定できるデータ、できないデータともに適用することができます。

クラスカル・ウォリス検定でデータ間に有意差が認められた場合、どの条件間に差があるのかを調べるために、下位検定として多重比較を行います。多重比較の方法は様々にありますが、わかりやすく適用範囲の広い Bonferroniの方法とその改良版、特に Shaffer の方法を使って p 値を調整しながらウィルコクソンの順位和検定を繰り返すのが良いと思います。

[2023/02/21 追記]
クラスカルウォリス検定を実施し、有意であれば多重比較を実施するという流れはメジャーだと思っていたのですが、始めから多重比較のみをすればよいという考えもあるようです。たしかに、一要因の場合には、多重比較でわかること ⊇ クラスカルウォリス検定でわかること、が多くの場合に成り立つような気がします。どちらがより適切かは、水準の数など実験デザインにも依存するように思うので、ケースバイケースで検討するしかないかなと考えています。

R Sample Code

# ライブラリのインポート
library(tidyverse)

# csvファイルの読み込み
# 整然データを想定: 一行が一観測。 "参加者ID, 独立変数, 従属変数" の形。
data <- read_csv("csvファイルのパス", show_col_types=F) %>%
    # コードの再利用可能性を高めるために列名を変更
    mutate(s=factor("参加者IDの列名"), A=factor("独立変数の列名"), y="従属変数の列名") %>%
    # 必要な列だけ抽出
    select(s, A, y) %>%
    # ソート
    arrange(A, s)

# 欠損値や外れ値など、解析から除外するべきデータがあれば削除
# 割愛

# 記述統計量の確認
data %>%
    group_by(A) %>%
    summarize(N=n(), Mean=mean(y), SD=sd(y), SE=sd(y)/sqrt(n()), Min=min(y), Median=median(y), Max=max(y))

# 簡単なグラフを描画して確認
data %>%
    # ggplot で作図開始。x軸は実験条件、y軸は評価指標（平均値）。
    ggplot(aes(x=A, y=y, fill=A)) +
    # 箱ひげ図
    geom_boxplot() +
    # 平均値
    stat_summary(fun=mean, geom="point", color="black", shape=4, size=2) +
    # 軸ラベルの設定
    labs(x="Condition", y="Measurement") +
    # カラーパレット
    scale_fill_brewer(palette="Blues") +
    # テーマの調整
    theme_light() +
    theme(legend.position="none") +
    theme(panel.grid.major.x=element_blank(), panel.grid.minor.y=element_blank())

# 等分散性の検定：p > alpha なら ok、そうでないならストップ。
library(car)
leveneTest(y~A, data=data)

# クラスカル・ウォリス検定の実行
(table <- kruskal.test(y~A, data=data))

# 効果量の計算
# 提示の意味はあまりなく、下位検定の効果量を提示すべきという説あり
# 案1 Morse(2009), Nye & Hans-Vaughn(2011) η2=χ2/(サンプルサイズ-1)
eta_squared <- table$statistic / (length(data$y)-1)
names(eta_squared) <- NULL
cat(paste("eta squared = ", eta_squared, "\n", sep=""))
# 案2 Field(2005, 2009) r=z/sqrt(N)
p <- table$p.value 
z <- qnorm(1 - p/2)#両側検定の場合
r <- z / sqrt(length(data$y))
cat(paste("r = ", r, "\n", sep=""))

# 下位検定：主効果が有意であった場合に実行
# Shaffer の方法で p 値を調整しながら、3回ウィルコクソンの順位和検定を行う例
library(exactRankTests)
a1 <- data %>% filter(A==levels(A)[1]) %>% pull(y)
a2 <- data %>% filter(A==levels(A)[2]) %>% pull(y)
a3 <- data %>% filter(A==levels(A)[3]) %>% pull(y)
t12 <- wilcox.exact(a1, a2, paired=F)
t13 <- wilcox.exact(a1, a3, paired=F)
t23 <- wilcox.exact(a2, a3, paired=F)
cohen.r <- function(x,y){
    t <- wilcox.exact(x,y,paired=F)
    p <- t$p.value
    z <- qnorm(1 - p/2)
    return (z / sqrt(length(c(x,y))))
    }
data.frame(term=c("a1.a2", "a1.a3", "a2.a3")) %>%
    mutate(v=c(t12$statistic, t13$statistic, t23$statistic)) %>%
    mutate(p=c(t12$p.value, t13$p.value, t23$p.value)) %>%
    mutate(r=c(cohen.r(a1,a2), cohen.r(a1,a3), cohen.r(a2,a3))) %>%
    # p 値の小さい順に並べる
    arrange(p) %>%
    # 各ステップで正味の仮説数はいくつか書きだす
    # ここでは主効果が有意という結果を使って、1ステップ目の正味の仮説数を減らしています
    mutate(k=c(1, 1, 1)) %>%
    # 各ステップの正味の仮説数を使って調整した p 値を判定に用いる（有意水準alphaの方を調整するという考え方もある）
    mutate(p.adj=p*k) %>%
    # p.adj < alpha ならば有意
    mutate(sig.=symnum(p.adj, corr=FALSE, na=FALSE, cutpoints=c(0,0.001,0.01,0.05,0.1,1), symbols=c("***","**","*","."," ")))

How to Report : LaTeX Sample

Japanese

等分散性 (F test, \textit{p} $>$ .05) が棄却されなかったため，クラスカル・ウォリス検定を実施した結果，yに対する要因Aの主効果は有意であった ($\chi^2$ (X) = X.XX, \textit{p} = .XX).
主効果が有意であったため、下位検定として Shaffer の方法のもとでウィルコクソンの順位和検定を繰り返し行った結果、yは条件A1が条件A2よりも有意に高かったが (\textit{p} = .XX, Cohen's \textit{r} = X.XX)、条件A1とA3間 (\textit{p} = .XX, Cohen's \textit{r} = X.XX) および条件A2とA3間 (\textit{p} = .XX, Cohen's \textit{r} = X.XX) には有意差は認められなかった。

English

Because homogeneity of variance assumption was not violated (F test, \textit{p} $>$ .05), we conducted a Kruskal-Wallis test. The result showed a significant main effect of A. Then, we conducted three Wilcoxon rank sum tests (Shaffer-corrected), which showed that Y was significantly higher in the A1 condition than in the A2 condition (\textit{p} = .XX, Cohen's \textit{r} = X.XX), whereas there were no significant differences between the A1 and A3 conditions (\textit{p} = .XX, Cohen's \textit{r} = X.XX) and between the A2 and A3 conditions (\textit{p} = .XX, Cohen's \textit{r} = X.XX).

For Your Information

効果量について

クラスカル・ウォリス検定の効果量については、定まっていないようです。いくつかの提案を見つけたので、サンプルコードに載せておきました。一方で、結局は主効果が有意であれば下位検定をすることになりますし、知ってうれしいのは下位検定で比較する条件間の差の大きさなので、下位検定の効果量を報告しておけばよいという考えもあるようです。

等分散性が仮定できない場合

等分散性が仮定出来ない場合には、クラスカル・ウォリス検定では精度が下がるそうです。これは、ウィルコクソンの順位和検定が、実は等分散性の仮定を必要とすることと同様です。

等分散性が仮定出来ない場合の検定方法として、一つはウェルチの一元配置分散分析を使うという手があるかもしれません。つまり、分散分析は正規性に関してロバストであるという前提で分散分析を強行するということです。

等分散性に関する脆弱性はあまり知られていないようなので、科学的に正しいかどうかを抜きにして考えれば、等分散性が仮定できないデータに対してクラスカル・ウォリス検定を適用してしまっても指摘されることはほとんどないかもしれません。根本的に解決したいと考えるなら、様々なデータの分布を仮定できる一般化線形混合効果モデルなどを勉強すれば、何か方法があるかもしれません。頑張りたいです。