分散分析と F 検定の関係が分からなくなることが多いのでメモ

二元配置分散分析（Two-Way ANOVA）

分散分析は、3 つ以上の群の平均値に差があるかどうかを評価する統計手法で、異なる条件や処理が平均値に与える影響を評価できる。なお、二元配置分散分析は、要因数が 2 種類の分散分析。

各々の群について、処理の影響を示す群間分散と処理と関係のない群内分散を求める。次に、群内分散に対する群間分散の比を求める。そして、この比を F 値として、F 分布から p 値を求める。このとき、二元配置分散分析は、異なる条件や処理が同時に作用したときの効果（交互作用）を評価することが特徴になる。なお、このときの帰無仮説は「すべての群の平均値が等しい」になる。

対応のない二元配置分散分析

一元配置分散分析は 1 要因の各水準からなる群に対して要因分散と誤差分散を求めて、誤差分散に対する要因分散の比が F 値になりました。

二元配置分散分析は、要因が 2 に増えますが、計算手順は一元配置分散分析と同じく誤差分散に対する要因分散の比を求めます。ただし、2 つの要因に加えて、その交互作用を計算する点が異なります。そして、$F_A$、$F_B$、$F_{AB}$のそれぞれについて F 検定を行い、有意差を確認します。

F 値の式

対応のない二元配置分散分析の場合、個体間の変動を考慮しない。

交互作用と誤差の自由度の計算に注意して、要因 A、要因 B、交互作用、誤差について、それぞれの変動、自由度から分散を計算できる。

項目	変動 ($SS$)	自由度 ($\nu$, $df$)	分散	F 値
全体	$$SS_T = \sum{i=1}^a \sum{j=1}^b \sum{k=1}^{n{ij}} (X{ijk} - \overline{X}{..})²$$	$$\nu_T = N - 1$$	$$-$$	$$-$$
要因 A	$$SS_A = \sum{i=1}^a n{i} (\overline{X}{i.} - \overline{X}{..})²$$	$$\nu_A = a - 1$$	$$MS_A = \frac{SS_A}{\nu_A}$$	$$F_A = \frac{MS_A}{MS_E}$$
要因 B	$$SS_B = \sum{j=1}^b n{j} (\overline{X}{.j} - \overline{X}{..})²$$	$$\nu_B = b - 1$$	$$MS_B = \frac{SS_B}{\nu_B}$$	$$F_B = \frac{MS_B}{MS_E}$$
交互作用	$$SS{AB} = \sum{i=1}^a \sum{j=1}^b n{ij} (\overline{X}{ij} - \overline{X}{i.} - \overline{X}{.j} + \overline{X}{..})²$$	$$\nu_{AB} = ab - 1 - \nu_A - \nu_B = \nu_A \times \nu_B$$	$$MS{AB} = \frac{SS{AB}}{\nu_{AB}}$$	$$F{AB} = \frac{MS{AB}}{MS_E}$$
誤差	$$SS_E = SS_T - SS_A -SS_B -SS{AB} \= \sum{i=1}^a \sum{j=1}^b \sum{k=1}^{n{ij}} (X{ijk} - \overline{X}_{ij})²$$	$$\nu_E = N - ab = \nu_T - \nu_A - \nu_B - \nu_{AB}$$	$$MS_E = \frac{SS_E}{\nu_E}$$	$$-$$

サンプル数
- $N$ $\cdots$ 全観測値の数
- $n_i$ $\cdots$ 要因 A の第 $i$ 水準に属する観測値の数
- $n_j$ $\cdots$ 要因 B の第 $j$ 水準に属する観測値の数
- $n_{ij}$ $\cdots$ 要因 A の第 $i$ 水準と要因 B の第 $j$ 水準に属する観測値の数 (繰り返し測定数)
水準数
- $a$ $\cdots$ 要因 A の水準数
- $b$ $\cdots$ 要因 B の水準数
観測値
- $X_{ijk}$ $\cdots$ 要因 A の第 $i$ 水準と要因 B の第 $j$ 水準に属する $k$ 番目の観測値
- $\overline{X}_{i.}$ $\cdots$ 要因 A の第 $i$ 水準に属する観測値の平均
- $\overline{X}_{.j}$ $\cdots$ 要因 B の第 $j$ 水準に属する観測値の平均
- $\overline{X}_{ij}$ $\cdots$ 要因 A の第 $i$ 水準と要因 B の第 $j$ 水準に属する観測値の平均
- $\overline{X}_{..}$ $\cdots$ 全観測値の平均

R で実装する

F 値と p 値は、対応無しの場合と同様に計算します。

デモデータの作成

デモ用のデータを作成します。このデータは、3 つの水準（A, B, C）を持つ要因 A（gA）と、2 つの水準（a, b）を持つ要因 B （gB）によって分類された測定値（value）からなります。

rm(list=ls(all=TRUE))
set.seed(123) # 乱数のシードを固定
data <- data.frame(
  gA = factor(rep(LETTERS[1:3], each = 4)), # 要因A、3水準
  gB = factor(rep(letters[1:2], each = 2)), # 要因B、2水準
  value = unlist(lapply(c(1, 4) * 40, sample, size = 6, replace = TRUE)) # 測定値
)

knitr::kable(sapply(levels(data$gA), function(i) {
  sapply(levels(data$gB), function(j) {
    paste(subset(data, subset = gA == i & gB == j)$value, collapse = ",")
  })
}))

	A	B	C
a	31,15	37,14	90,91
b	14,3	118,153	91,92

表の形は対応のある一元配置分散分析と同じ。

でも、対応のある一元配置分散分析が a と b は互いに影響しない(独立になる)ことに対して、対応のない二元配置分散分析は a と b は互いに影響しあう可能性、交互作用を考えることになる。

ここが大きく違う点で、要因A、要因Bに加えて、交互作用に対してもF検定をすることになる。

分散分析の関数を使わない

# 各水準の観測個数
nij <- table(data[,-3])

# 総変動
# 平均
mean_T <- mean(data$value)
# 変動、偏差平方和
ss_T <- sum((data$value - mean_T)^2)
# 自由度、総データ数 - 1
df_T <- nrow(data) - 1

# 要因A
# 平均
mean_A <- as.vector(tapply(data$value, data$gA, mean))
# 変動、偏差平方和
ss_A <- sum((mean_A[data$gA] - mean_T)^2)
# 自由度、水準数 - 1
df_A <- length(levels(data$gA)) - 1
# 分散
ms_A <- ss_A / df_A

# 要因B
# 平均
mean_B <- as.vector(tapply(data$value, data$gB, mean))
# 変動、偏差平方和
ss_B <- sum((mean_B[data$gB] - mean_T)^2)
# 自由度、水準数 - 1
df_B <- length(levels(data$gB)) - 1
# 分散
ms_B <- ss_B / df_B

# 交互作用
# 平均
mean_AB <- tapply(data$value, list(data$gA, data$gB), mean)
# 変動、偏差平方和
ss_AB <- sum(nij*t(t(mean_AB - mean_A) - mean_B + mean_T)^2)
# 自由度、(要因A水準数 - 1)(要因B水準数 - 1)
df_AB <- df_A*df_B
# 分散
ms_AB <- ss_AB / df_AB

# 誤差
# 変動、偏差平方和
ss_E <- ss_T - ss_A - ss_B - ss_AB
# ss_E <- sum((data$value - mapply(rep, t(mean_AB), t(nij)))^2)
# 自由度、(要因A水準数 - 1)(要因B水準数 - 1)
df_E <- df_T - df_A - df_B - df_AB
# 分散
ms_E <- ss_E / df_E

# F値、要因分散/誤差分散
f_A <- ms_A / ms_E
f_B <- ms_B / ms_E
f_AB <- ms_AB / ms_E

# F値からp値を求める
p_A <- pf(f_A, df_A, df_E, lower.tail = FALSE)
p_B <- pf(f_B, df_B, df_E, lower.tail = FALSE)
p_AB <- pf(f_AB, df_AB, df_E, lower.tail = FALSE)

# 結果の表示
cat("factor A - F-value:", round(f_A, 2), "/ p-value:", round(p_A, 4), "\n")
cat("factor B - F-value:", round(f_B, 2), "/ p-value:", round(p_B, 4), "\n")
cat("factor AB - F-value:", round(f_AB, 2), "/ p-value:", round(p_AB, 4), "\n")

factor A - F-value: 37.38 / p-value: 4e-04
factor B - F-value: 17.46 / p-value: 0.0058
factor AB - F-value: 25.9 / p-value: 0.0011

分散分析の関数を使う

二元配置分散分析は 2 要因を$\ast$で繋ぐ。

# 対応のある一元配置分散分析
res <- aov(value ~ gA * gB, data = data)
summary(res)

            Df Sum Sq Mean Sq F value  Pr(>F)
gA           2  13287    6644   37.38 0.00041 ***
gB           1   3104    3104   17.46 0.00582 **
gA:gB        2   9207    4604   25.90 0.00112 **
Residuals    6   1067     178
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

要因	F 値	p 値	$\alpha$ = 0.05	$\alpha$ = 0.01
要因 A	37.38	0.00041	有意差あり	有意差あり
要因 B	17.46	0.00582	有意差あり	有意差あり
交互作用	25.90	0.00112	有意差あり	有意差あり

つまり、いずれの要因も有意差があって、交互作用もある、という結果になる。

ただし、交互作用に有意差があっても、要因 A と要因 B のどの水準の組み合わせに有意差があるのかは分からない。

この先の解析は多重比較になる。

日々のつれづれ

不惑をむかえ戸惑いを隠せない男性の独り言

分散分析のメモ・その３