🎳

HCI研究者の心理統計ノート：一元配置分散分析 / One-way ANOVA

2023/01/09に公開

statistics

research

tech

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

Introduction

一元配置分散分析 (One-way ANOVA) は、3つ以上の対応のないデータ間の平均値に差があるかどうかを検定する方法で、データに正規性および等分散性が仮定できる場合に使用できます。対応のあるデータの場合には、繰り返しのある一元配置分散分析を参照してください。また、正規性が仮定できない場合は、クラスカル・ウォリス検定を参照してください。等分散性が仮定出来ない場合については、ウェルチの一元配置分散分析を参照してください。

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

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

分散分析を R で実行する方法はいくつもあります。（一元配置分散分析の場合は）どれを使っても問題ないと思いますし、いくつか使ってみて同じ結果になることを確認してみてもいいと思います。特に、井関龍太先生が公開してくださっているANOVA君 (anovakun) は、二元配置分散分析やさまざまな効果量の計算などにも対応しており、簡単に分散分析が実施できるので、おすすめです。しかし、はじめからすべてがブラックボックスになってしまうことに私は多少なりとも抵抗があるので、この記事のメインのコードでは anovakun 以外の実装方法を紹介します。anovakun を使った方法も、記事下部の For Your Information に紹介しています。

[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 %>%
    # 平均と標準誤差を計算
    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)

# 等分散性の検定 p > alpha なら ok (棄却された場合は、For Your Informationを参照)
library(car)
leveneTest(y~A, data=data)

# 一元配置分散分析
library(car)
options(contrasts=c("contr.sum", "contr.sum"))
(table <- Anova(lm(y~A, data=data), type=3))

# 効果量 partial eta squared
SS.effect.A <- table["A", "Sum Sq"]
SS.residuals <- table["Residuals", "Sum Sq"]
peta <- SS.effect.A / (SS.effect.A + SS.residuals)
cat(paste("partial eta squared = ", peta, "\n", sep=""))

# 下位検定：主効果が有意であった場合に実行
# Shaffer の方法で p 値を調整しながら、3回スチューデントの 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)
t12 <- t.test(a1, a2, var.equal=T)
t13 <- t.test(a1, a3, var.equal=T)
t23 <- t.test(a2, a3, var.equal=T)
data.frame(term=c("a1.a2", "a1.a3", "a2.a3")) %>%
    mutate(t=c(t12$statistic, t13$statistic, t23$statistic)) %>%
    mutate(df=c(t12$parameter, t13$parameter, t23$parameter)) %>%
    mutate(p=c(t12$p.value, t13$p.value, t23$p.value)) %>%
    # 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

正規性 (Shapiro-Wilk test, \textit{p} $>$ .05) および等分散性 (Levene 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), 条件A2とA3には有意差は認められなかった (\textit{t} (X) = X.XX, \textit{p} = .XX).

English

Because normality (Shapiro-Wilk's normality test, \textit{p} $>$ .05) and homogeneity of variance (Levene test, \textit{p} $>$ .05) assumptions were not violated, we conducted a 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 three Student'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 A2 and A3 conditions (\textit{t} (X) = X.XX, \textit{p} = .XX).

For Your Information

anovakun を使う方法とその他の実装方法

anovakun の使い方は、公式ページにとても丁寧に書いてくださっているので、詳しくはそちらを参照してください。

以下は、上記のサンプルコードとほとんど同じ結果を出力するサンプルコードです。ただし、分散分析の結果については一致するのですが、多重比較の数値が微妙に合わない場合があり困っています。anovakun を解読してみると、どうやら t 値の求め方が異なっているように思います。何が正しいのかよくわかっていませんが、メインのサンプルコードで示した方法の方がどんな下位検定をしたのか説明しやすいという側面はあるかもしれません。無責任で申し訳ないです。

source("~/Documents/R/anovakun_486.txt")
anovakun(data %>% select(A, y), "As", 3, peta=T, fs1=T)

それから、インターネット上で見かけたその他の実装方法についても、せっかくなので記しておきます。一元配置分散分析では、すべて同じ結果になると思います。ただし、二元配置分散分析以上になると、主に平方和の計算方法の違いに起因して、適切かどうかがわかれてきます。サンプルコードに示した方法は、平方和を明示的に指定できるものにしました。

anova(lm(y~A, data=data))
anova(aov(y~A, data=data))
oneway.test(y~A, data=data, var.equal=T)
summary(aov(y~A, data=data))

効果量について

分散分析の効果量は、以下の式で定義される偏イータ二乗を使っておくのがよいと思われます。

{\eta_p}^2 = \frac{（ある効果の平方和）}{（ある効果の平方和 + 対応する誤差平方和）} = \frac{SS_{effect}}{SS_{effect} + SS_{residuals}}

これは、一元配置分散分析の場合には、イータ二乗と同じ値になります。ちなみに、イータ二条の定義は以下の通りです。

\eta^2 = \frac{（ある効果の平方和）}{（全体の平方和）} = \frac{SS_{effect}}{SS_{Total}}

他にも、一般化イータ二乗、オメガ二乗、などなどが提案されているようですが、よく見るのは偏イータ二乗です。ただし、偏イータ二乗には大きさの目安が提案されていないため、その点についてはやや扱いにくいかもしれません。一元配置分散分析の場合は、偏イータ二乗とイータ二乗が同じなので、イータ二乗の目安に倣って考えることが出来ます (小: 0.01, 中: 0.06, 大: 0.14)。しかし、基本的には、そうした目安に基づいた議論よりも、同分野の先行研究の効果量と比較した議論のほうが有意義だと思います。

多重比較の方法について

Bonferroni の方法と、その改良版である Holm, Shaffer の方法などについては、以下のANOVA君のページがとても参考になるので、ぜひ一読することをおすすめします。他にも Tukey, Dunnett, Williamsなどの手法がありますが、Bonferroni 系とどのような使い分けをすればよいのかについては、私は今ひとつわかっていません。一つは、Bonferroni 系ならば、すべての組み合わせを比較したり、統制群とその他すべての条件を比較したり、といった決まった形に縛られず、行う検定の回数のみから p 値を調整できるという利点があります。個人的には、ひとまずは、適用可能な場合が多く、わかりやすい Bonferroni 系、特に検出力を高めた Shaffer の方法を使っておけばよいと考えています。

Shafferの方法において有意水準を割る（または p 値にかける）正味の仮説数については、Holland & Copenhaver (1987)に表が掲載されています。3群の比較であれば k = (3, 1, 1), 4群の比較であれば k = (6, 3, 3, 3, 2, 1) で、p.adj = p * k[i] を有意水準 alpha と比較する、または p を有意水準 alpha / k[i] と比較して判定することになります。分散分析の主効果が有意であれば、1ステップ目の正味の仮説数を2ステップ目と同じまで減らすことが出来るため、3群の比較であれば k = (1, 1, 1), 4群であれば k = (3, 3, 3, 3, 2, 1) で良いということになります。Bonferroni 系は、比較する回数が多くなるほど検出力が弱まるため、5群、6群とある場合には、他の多重比較方法の使用を検討するのが良いかもしれません。