🍂

HCI研究者の心理統計ノート：繰り返しのある一元配置分散分析 / One-way repeated measures ANOVA

2023/01/10に公開

statistics

research

tech

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

Introduction

繰り返しのある一元配置分散分析 (One-way repeated measures ANOVA) は、3群以上の対応のあるデータ間の平均値に差があるかどうかを検定する方法で、データに正規性および球面性が仮定できる場合に使用できます。対応のないデータの場合には一元配置分散分析、正規性が仮定できない場合はフリードマン検定を参照してください。球面性の仮定については、球面性の検定を行い、満たされなければ補正を行います。

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

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

分散分析を R で実行する方法はいくつかあり、特に井関龍太先生が公開してくださっているANOVA君 (anovakun) は、二元配置分散分析やさまざまな効果量の計算などにも対応しており、簡単に分散分析が実施できるのでおすすめです。しかし、はじめからブラックボックスにしすぎてしまうことに私は多少なりとも抵抗があり、この記事のメインのコードでは anovakun 以外の実装方法を紹介します。anovakun を使った方法は、記事下部の For Your Information に紹介しています。私はこれらの方法でダブルチェックすることが多いです。慣れてくれば、というか、本当は最初から anovakun の恩恵に預からせてもらっておけばOKだと思います。

[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)

# 球面性の検定と繰り返しのある一元配置分散分析
# Mauchlyの球面性検定が有意でなければOK。
# 有意なら、球面性の補正の出力部分を確認して、Geenhouse-Geisser(GG)またはHuynh-Feldt(HF)の補正のイプシロンの値と補正したp値を確認する。
library(ez)
options(contrasts=c("contr.sum", "contr.poly"))
(table <- ezANOVA(data=data, dv=.(y), wid=.(s), within=.(A), type=3))

# 効果量 partial eta squared
F <- table$ANOVA[["F"]]
Df.A <- table$ANOVA[["DFn"]]
Df.R <- table$ANOVA[["DFd"]]
peta <- 1 / (1 + Df.R/(F*Df.A))
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, paired=T)
t13 <- t.test(a1, a3, paired=T)
t23 <- t.test(a2, a3, paired=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) および球面性 (Mauchly 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 sphericity (Mauchly 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 paired \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 の使い方は、公式ページにとても丁寧に書いてくださっているので、詳しくはそちらを参照してください。

以下は、上記のサンプルコードとほとんど同じ結果を出力するサンプルコードです。ただし、球面性の検定として Mauchly の検定ではなく Mendoza の多標本球面性検定を実施します。オプションを指定すれば、Mauchly の検定にも変更出来ます。

source("~/Documents/R/anovakun_486.txt")
anovakun(data %>% pivot_wider(names_from="A", values_from="y") %>% select(-s), "sA", length(levels(data$A)), peta=T, fs1=T)

球面性の仮定について

球面性の仮定とはなにか、球面性の検定について、補正とは、などなど、井関先生がANOVA君のサイトにとてもわかりやすく丁寧に書いてくださっています。

同じく、Mauchly の検定と Mendoza の検定の違いについても書かれています。これを読む限り、anovakun のデフォルトがそうであるように、Mendoza の検定を使っていた方が適切なような気もしてきます。

その他

効果量や多重比較などについて、共通する情報があるので、一元配置分散分析のページも確認してみてください。