Zenn
🍱

HCI研究者の心理統計ノート：Two-way ANOVA for between-subject design

2023/05/29に公開
R
statistics
research
tech

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

https://zenn.dev/tmizuho/articles/e8890371f683bd

Introduction

参加者間計画における二元配置分散分析 (Two-way ANOVA for between-subject design) は、参加者間要因 × 参加者間要因の二要因参加者間計画において、各要因の主効果と要因間の交互作用が有意であるかどうかを検定する方法です。正規性および等分散性が仮定出来る場合に使用できます。

参加者間要因 × 参加者内要因の場合は 混合計画における二元配置分散分析 (Two-way ANOVA for mixed design)、参加者内要因 × 参加者内要因の場合は 参加者内計画における二元配置分散分析 (Two-way ANOVA for within-subject design) を参照してください。正規性が仮定出来ない場合には、参加者間計画における整列ランク変換を施した二元配置分散分析 (Two-way ANOVA with ART for between-subject design) を参照してください。

HCI研究では、要因A（水準a1, a2, ..., an） × 要因B（水準b1, b2, ..., bm）の n*m 個のグループに参加者を分けて、従属変数を比較する実験で使用します。一元配置分散分析と大きく違う点は、要因間の交互作用を知ることが出来る点です。例えば、手法（既存 vs 提案） × 性別（男性 vs 女性）の実験計画では、男性参加者には既存手法の方が有効だが女性参加者には提案手法の方が有効である、などといった２つの要因が複合的に評価指標に作用する可能性を検証することができます。

分散分析を R で実行する方法はいくつもありますが、二元配置分散分析の場合には平方和の計算方法によって適切な方法とそうでない方法があるので注意が必要です。非釣り合い型データ（グループごとの人数が異なる）の場合に、平方和の計算方法に依存して結果が変わってしまう可能性があります。詳細は井関龍太先生のページがとてもわかりやすいので確認してください。結論としては、タイプIII平方和を使っておけば安心だと思っておけばよいと思われます。

最も簡単な実行方法は、井関先生が公開してくださっている ANOVA君 (anovakun) を利用する方法です。しかし、はじめからすべてがブラックボックスになってしまうことに私は多少なりとも抵抗があり、この記事のメインのサンプルコードでは anovakun 以外の実装方法を紹介します。anovakun を用いる方法も記事下部の For Your Information に紹介しています。私はこれらを用いてダブルチェックしています。

二元配置分散分析から始まる一連の統計フローは次の通りです。

  1. 要因A×要因Bの二元配置分散分析を実施する
  2. 交互作用が有意であれば、単純主効果の検定を実施する（e.g., A=aiのグループのみで、要因Bについての一元配置分散分析または t 検定を行う。B=bjのグループのみで、要因Aについての一元配置分散分析または t 検定を行う。）
    a. 単純主効果が有意であれば、p 値を補正しながら多重比較を実施して、終わり
    b. 単純主効果が有意でなければ、終わり
  3. 交互作用が有意でなく要因Aの主効果が有意であれば、要因Aのどの水準間に差があるのか調べるために、p 値を補正しながら多重比較を実施する
  4. 交互作用が有意でなく要因Bの主効果が有意であれば、要因Bのどの水準間に差があるのか調べるために、p 値を補正しながら多重比較を実施する
  5. 交互作用が有意でなく、いずれの主効果も有意でなければ、終わり

anovakun は自動でこのフローを実行してくれるので、とてもありがたいです。しかし、こういうフローになっているということは理解しておく必要があります。

特に、交互作用が有意である場合は、主効果が有意であってもあまり意味がありません。この点を勘違いしている学生も割といるような気がします。交互作用が有意である場合は、単純主効果の検定を行うことによって、データの傾向を正しく把握し論じるようにしましょう。

以下のサンプルコードでは、2×3の2要因参加者間計画において、交互作用が有意で単純主効果も有意であるパターンを載せておきます。このパターンを抑えておけば、その他の実験デザインの場合にも拡張しやすいのではないかと思います。とはいえ、anovakun を使えば、すべてのパターンを自動で判別して実行してくれます。

R Sample Code

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

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

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

# 記述統計量の確認
data %>%
    group_by(A, B) %>%
    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, B) %>%
    summarize(Mean=mean(y), SE=sd(y)/sqrt(n())) %>%
    # ggplot で作図開始。x軸は要因A、y軸は評価指標（平均値）、凡例は要因B。
    ggplot(aes(x=A, y=Mean, fill=B)) +
    # 棒グラフ
    geom_bar(stat="identity", width=0.9, position=position_dodge(0.9)) +
    # エラーバー
    geom_errorbar(aes(ymin=Mean-SE, ymax=Mean+SE), width=.2, size=.75, position=position_dodge(0.9)) +
    # 軸ラベルの設定
    labs(x="Factor A", y="Measurement", fill="Factor B") +
    # カラーパレット
    scale_fill_brewer(palette="Blues") +
    # テーマの調整
    theme_light() +
    theme(legend.position="top") +
    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*B, data=data)

# ANOVA 
library(car)
options(contrasts=c("contr.sum", "contr.poly"))
(table <- Anova(lm(y~A*B, data=data), type=3))

# 効果量 partial eta squared
SS.effect.A <- table["A", "Sum Sq"]
SS.effect.B <- table["B", "Sum Sq"]
SS.effect.AxB <- table["A:B", "Sum Sq"]
SS.residuals <- table["Residuals", "Sum Sq"]
peta.A <- SS.effect.A / (SS.effect.A + SS.residuals)
peta.B <- SS.effect.B / (SS.effect.B + SS.residuals)
peta.AxB <- SS.effect.AxB / (SS.effect.AxB + SS.residuals)
data.frame(peta.A=peta.A, peta.B=peta.B, peta.AxB=peta.AxB)

# ここまでは水準数に依らず、適用可能
# ここからは A (3水準; a1, a2, a3) × B (2水準; b1, b2) の実験デザインを想定したサンプルコード
# 交互作用が有意であれば、単純主効果の検定を行う
# (1-1) B=b1のデータのみを抽出して、要因Aについて一元配置分散分析を実施する
(table <- Anova(lm(y~A, data=data %>% filter(B==levels(B)[1])), type=3))
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=""))

# (1-2) B=b1のもとでの要因Aの単純主効果が有意であれば、Shaffer の方法で p 値を調整しながら、3回スチューデントの t　検定を行う
# ウェルチの t 検定の方がいいかもしれない。その場合は var.equal=F に書き換える。
a1 <- data %>% filter(B==levels(B)[1], A==levels(A)[1]) %>% pull(y)
a2 <- data %>% filter(B==levels(B)[1], A==levels(A)[2]) %>% pull(y)
a3 <- data %>% filter(B==levels(B)[1], 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)) %>%
    # 各ステップの正味の仮説数 k を使って p 値 を調整する p_adj=p*k 
    # または、alphaの方を alpha_adj=alpha/k と調整してもよい
    mutate(p.adj=p*k) %>%
    # p.adj < alpha （または、p < alpha_adj） なら有意
    mutate(sig.=symnum(p.adj, corr=FALSE, na=FALSE, cutpoints=c(0,0.001,0.01,0.05,0.1,1), symbols=c("***","**","*","."," ")))
    
# (2) B=b2のデータについても、同様に要因Aの単純主効果の検定を行う。必要に応じて、下位検定を行う。
# 割愛

# (3) A=a1のデータのみを抽出して、要因Bの単純主効果の検定を行う
# 要因Bは二水準なので、t検定を実施する
b1 <- data %>% filter(A==levels(A)[1], B==levels(B)[1]) %>% pull(y)
b2 <- data %>% filter(A==levels(A)[1], B==levels(B)[2]) %>% pull(y)
(t12 <- t.test(b1, b2, var.equal=T)

# (4) A=a2のデータについても同様に、要因Bの単純主効果の検定を行う。（割愛）
# (5) A=a3のデータについても同様に、要因Bの単純主効果の検定を行う。（割愛）

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), 〇〇の主効果 (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX), および交互作用 (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX) が有意であった．
交互作用が有意であったため，単純主効果の検定を実施した．一元配置分散分析の結果，条件B1では要因Aの単純主効果が有意であった (\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).他方，条件B2では要因Aの単純主効果が有意ではなかった (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX).
また，要因Bの単純主効果の検定としてスチューデントの \textit{t} 検定を実施した結果，条件A1, A2, A3 のいずれでも要因Bの単純主効果は有意でなかった (\textit{p}s > .1).

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 two-way ANOVA. The result showed a significant main effect of 〇〇 (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX). The main effect of 〇〇 was also significant (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX). Furthermore, the interaction effect of A and B was significant (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX).
Then, we conducted two one-way ANOVA to verify the simple main effects of A. In the B1 condition, the result showed a significant simple-main effect of A (\textit{F} (X, XX) = X.XX, \textit{p} = .XX, ${\eta_p}^2$ = 0.XX), while no significant simple-main effect was observed in the B2 condition. Three Student's \textit{t} tests (Shaffer-corrected) showed that y in the A1-B1 condition was significantly higher than in the A2-B1 condition (\textit{t} (X) = X.XX, \textit{p} = .XX) and in the A3-B1 conditon (\textit{t} (X) = X.XX, \textit{p} = .XX), but there was no significant difference between the A2-B1 and A3-B1 conditions (\textit{t} (X) = X.XX, \textit{p} = .XX). 
Additionally, we performed three Student's \textit{t} tests to examine the simple main effects of B, which showed no significant differences between the A1-B1 and A1-B2 conditions, between the A2-B1 and A2-B2 conditions, nor between the A3-B1 and A3-B2 conditions (\textit{p}s $>$ .1).

For Your Information

anovakun を使う方法

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

anovakun を使わない実装方法をメインに示しましたが、正直 anovakun を使ったほうが簡単です。単純主効果の検定や多重比較について、特に恩恵を感じます。本当にありがとうございます。まあ、面倒でなければ、いろいろな方法でダブルチェックするのは良いことだと思います。

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

source("~/Documents/R/anovakun_486.txt")
anovakun(data %>% select(A, B, y), "ABs", length(levels(data$A)), length(levels(data$B)), peta=T, fs1=T)

効果量について

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

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

他にも、一般化イータ二乗、オメガ二乗、などなどが提案されているようですが、よく見るのは偏イータ二乗です。ただし、偏イータ二乗には大きさの目安が提案されていないため、その点についてはやや扱いにくいかもしれません。個人的には、イータ二乗と同じ目安(小: 0.01, 中: 0.06, 大: 0.14)を使って「効果量は大きめかもしれないな」ぐらいは考えることがありますが、論文には書けません。

等分散性の仮定について

等分散性の仮定について、釣り合い型ならばロバスト、とどこかで見たような気がしますが、まだちゃんと調査していません。等分散性が棄却された場合に、適用できる代替手法の検討がついていないため、すごく困ります。一般化線形モデルなどの上位の概念を勉強する必要がありそうです。

Discussion