Zenn
🎍

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

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

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

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

Introduction

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

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

HCI研究では、参加者間要因A（水準a1, a2, ..., an） × 参加者内要因B（水準b1, b2, ..., bm）の n*m 個の条件で、従属変数を比較する実験で使用します。例えば、手法（既存 vs 提案） × 計測タイミング（直後 vs 1週間後）の実験計画では、直後は条件間に差がないが1週間後に計測すると提案手法の方が有効である、などといった２つの要因が複合的に評価指標に作用する可能性を検証することができます。

分散分析を 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 は自動でこのフローを実行してくれるので、とてもありがたいです。しかし、こういうフローになっているということは理解しておく必要があります。

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

以下のサンプルコードでは、分散分析を実施するところまでをカバーしています。下位検定として、単純主効果の検定や多重比較を行う部分については、水準数や主効果・交互作用の有意・非有意に依存して数多のパターンが存在するので、割愛しています。参加者間計画における二元配置分散分析 (Two-way ANOVA for between-subject design) の記事では、割愛せずに具体例を記述しようと努めているので、下位検定のイメージがつかめない人は参照してください。

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 1", y="Measurement", fill="Factor 2") +
    # カラーパレット
    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
library(car)
leveneTest(y~A, data=data)#Bについては、参加者内要因なので、検定しなくていいはず

# 球面性の検定 + ANOVA + 効果量 partial eta squared
library(ez)
options(contrasts=c("contr.sum", "contr.poly"))
table <- ezANOVA(data=data, dv=.(y), wid=.(s), between=.(A), within=.(B), type=3)
table$ANOVA$peta <- with(table$ANOVA, 1 / (1 + DFd/(F*DFn)))
table

# 交互作用が有意であれば、単純主効果の検定 -> 多重比較 を行う
# (A) 参加者間要因Aについては、一元配置分散分析やウェルチの t 検定
# (B) 参加者内要因Bについては、繰り返しのある一元配置分散分析や対応のある t 検定
# を行うことになる。割愛。

# 交互作用が有意でなく、主効果が有意であれば、その要因について多重比較を行う
# 割愛

How to Report : LaTeX Sample

Japanese

正規性 (Shapiro-Wilk test, \textit{p} $>$ .05) および等分散性 (Levene test, \textit{p} $>$ .05), 球面性 (Mauchly 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) が有意であった．
交互作用が有意であったため，単純主効果の検定を行った．（以下略）

English

Because normality (Shapiro-Wilk's normality test, \textit{p} $>$ .05), homogeneity of variance (Levene test, \textit{p} $>$ .05), and sphericity (Mauchly 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 ...

For Your Information

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

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

メインのサンプルコードでは anovakun を使わない実装方法を示しましたが、正直 anovakun を使ったほうが簡単です。単純主効果の検定や多重比較について、特に恩恵を感じます。本当にありがとうございます。以下は、上記のサンプルコードとほとんど同じ結果を出力するサンプルコードです。ただし、球面性の検定として Mauchly の検定ではなく Mendoza の多標本球面性検定を実施します。オプションを指定すれば、Mauchly の検定にも出来ます。

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

球面性の仮定について

球面性の仮定とはなにか、球面性の検定について、補正とは、などなど、井関先生がANOVA君のサイトにとてもわかりやすく丁寧に書いてくださっています。
http://riseki.php.xdomain.jp/index.php?ANOVA君/球面性検定の出力

同じく、Mauchly の検定と Mendoza の検定の違いについても書かれています。これを読む限り、anovakun のデフォルトがそうであるように、Mendoza の検定を使っていた方が適切なような気もしてきます。
http://riseki.php.xdomain.jp/index.php?ANOVA君/多標本球面性への対応

その他

共通する情報があるので、参加者間計画の二元配置分散分析のページも確認してみてください。
https://zenn.dev/tmizuho/articles/4f33cb07840511

Discussion