🍔

HCI研究者の心理統計ノート：ウェルチの一元配置分散分析 / Welch one-way ANOVA

2023/01/09に公開

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

Introduction

ウェルチの一元配置分散分析 (Welch one-way ANOVA) は、3つ以上の対応のないデータ間の平均値に差があるかどうかを検定する方法で、データに正規性が仮定できる場合に使用できます。通常の一元配置分散分析と比べて、等分散性の仮定が必要ない点が異なります。等分散性が仮定できるデータ、できないデータともに適用することができます。

HCI研究では、1要因参加者間計画 (between-subject design) の実験において、条件 1,2,...k 間で評価指標に統計的に有意な差があるかどうかを調べる際に使用します (e.g., 既存手法1 vs 既存手法2 vs 提案手法)。その要因が評価指標に与える効果を主効果 (main effect) と言い、一元配置分散分析によって帰無仮説が棄却された場合には"主効果が有意"、"有意な主効果が認められた"などと表現します。

分散分析で主効果が有意であった場合、どの条件の組み合わせに差があるのかを調べるために下位検定を実施します。この下位検定は多重比較と呼ばれ、様々な方法がありますが、適用範囲が広く使いやすいのは Bonferroni の方法とその改良版を用いて p 値を調整しながら t 検定を繰り返す方法だと思います。

一元配置分散分析を R で実行する方法はいくつもありますが、ウェルチの方法を実行できる方法は限られているように思います。ここでは one-way.test() 関数でオプション引数を指定する方法を紹介します。ウェルチの分散分析は、そもそもあまりメジャーではないのかもしれません。しかし、分散分析は正規性については頑健であっても、（特に釣り合い型でない場合は）等分散性については頑健でないようなので (Zimmer, 2004)、等分散性が仮定できない場合にはウェルチの方法を使うべきだと思われます。さらに言えば、分散分析の前に等分散性の検定を行うことは検定の多重性の問題を生じるとの考えから、等分散性の検定を行わずにウェルチの一元配置分散分析を実施した方がよい、との提案もあるようです (井口, 2018)。

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

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 %>%
    # 平均と標準誤差を計算
    group_by(A) %>%
    summarize(Mean=mean(y), SE=sd(y)/sqrt(n())) %>%
    # ggplot で作図開始。x軸は実験条件、y軸は評価指標（平均値）。
    ggplot(aes(x=A, y=Mean, fill=A)) +
    # 棒グラフ
    geom_bar(stat="identity") +
    # エラーバー
    geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2, size=.75) +
    # 軸ラベルの設定
    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 (棄却された場合は、クラスカル・ウォリス検定へ)
shapiro.test(data$y)

# ウェルチの一元配置分散分析
(table <- oneway.test(y~A, data=data, var.equal=F))

# 効果量 partial eta squared
F <- table$statistic[["F"]]
Df.A <- table$parameter[["num df"]]
Df.R <- table$parameter[["denom df"]]
peta <- 1 / (1 + Df.R/(F*Df.A))
cat(paste("partial eta squared = ", peta, "\n", sep=""))

# 下位検定：主効果が有意であった場合に実行
# Shaffer の方法で p 値を調整しながら、4群間で（6回の）ウェルチの t　検定を行う例
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)
a4 <- data %>% filter(A==levels(A)[4]) %>% pull(y)
t12 <- t.test(a1, a2, var.equal=F)
t13 <- t.test(a1, a3, var.equal=F)
t14 <- t.test(a1, a4, var.equal=F)
t23 <- t.test(a2, a3, var.equal=F)
t24 <- t.test(a2, a4, var.equal=F)
t34 <- t.test(a3, a4, var.equal=F)
data.frame(term=c("a1.a2", "a1.a3", "a1.a4", "a2.a3", "a2.a4", "a3.a4")) %>%
    mutate(t=c(t12$statistic, t13$statistic, t14$statistic, t23$statistic, t24$statistic, t34$statistic)) %>%
    mutate(df=c(t12$parameter, t13$parameter, t14$parameter, t23$parameter, t24$parameter, t34$parameter)) %>%
    mutate(p=c(t12$p.value, t13$p.value, t14$p.value, t23$p.value, t24$p.value, t34$p.value)) %>%
    # p 値の小さい順に並べる
    arrange(p) %>%
    # 各ステップで正味の仮説数はいくつか書く
    # ここでは主効果が有意という結果を使って正味の仮説数をさらに減らす方法を使っています
    mutate(k=c(3, 3, 3, 3, 2, 1)) %>%
    # 各ステップの正味の仮説数を使って p 値を調整する（有意水準の方を調整するという考え方もある）
    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

正規性が棄却されなかったため (Shapiro-Wilk test, \textit{p} $>$ .05)，ウェルチの一元配置分散分析を実施した結果，〇〇の主効果は有意であった (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX). 
下位検定として Shaffer の方法で p 値を調整しながらウェルチの t 検定を繰り返し行った結果，条件A1は条件A2よりも有意に高く (\textit{t} (X) = X.XX, \textit{p} = .XX), 条件A1は条件A3よりも有意に高かったが (\textit{t} (X) = X.XX, \textit{p} = .XX), それ以外の組み合わせには有意差は認められなかった (\textit{p}s $>$ .1).

English

Because normality assumption was not violated (Shapiro-Wilk's normality test, \textit{p} $>$ .05), we conducted an Welch's one-way ANOVA. The result showed a significant main effect of 〇〇 (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX).
Then, we conducted six Welch's \textit{t} tests (Shaffer-corrected), which showed that y in the A1 condition was significantly higher than in the A2 condition (\textit{t} (X) = X.XX, \textit{p} = .XX) and in the A3 conditon (\textit{t} (X) = X.XX, \textit{p} = .XX), but there was no significant difference between the other conditions (\textit{p}s $>$ .1).

For Your Information

一元配置分散分析のページも参照してください

効果量や多重比較などについては、多くの説明が共通するので、以下のページを参照してください。

正規性が棄却された場合は？

一元配置分散分析のノンパラメトリックバージョンは、クラスカル・ウォリス検定 (Kruskal-Walis test) です。しかし、クラスカル・ウォリス検定も等分散性が仮定できない場合に検定の精度が低下するようです。したがって、正規性が仮定できず、等分散性も仮定できない場合には、分散分析が正規性に関してロバストであると主張して、ウェルチの分散分析をするというのも一つの案になるかもしれません。

また、そうした危うい橋を渡るよりは、分散分析というしがらみから出て、任意の事前分布を仮定できる一般化線形混合効果モデルなどの上位互換を一生懸命勉強するのが良いような気もしています。いつか使いこなせるようになったら記事を書きたいと思います。