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

一元配置分散分析

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

各々の群について、処理の影響を示す群間分散と処理と関係のない群内分散を求める。次に、群内分散に対する群間分散の比を求める。そして、この比を F 値として、F 分布から p 値を求める。なお、このときの帰無仮説は「すべての群の平均値が等しい」になる。

対応がない場合の一元配置分散分析

同じ個体や対象から複数回測定されたものではないとき、「対応がない」になる。

F 値の式

\begin{eqnarray}

F_{(\nu_1,\nu_2)} &=& \frac{要因分散}{誤差分散}
&=& \frac{\frac{群間変動}{群間変動の自由度}}{\frac{群内変動}{群内変動の自由度}}
&=& \frac{\frac{反復数 \times \sum(群平均 - 総平均)^2}{処理数 - 1}}{\frac{\sum\sum(測定値 - 群平均)^2}{データ総数 - 処理数}}
&=& \frac{\frac{n\sum(\bar{x_i}-\bar{\bar{x}})^2}{\nu_1}}{\frac{\sum\sum(x-\bar{x_i})^2}{\nu_2}}

\end{eqnarray}

$F_{(\nu_1,\nu_2)}$F は、群間の平均値のばらつき（要因分散）が、群内のデータのばらつき（誤差分散）に比べてどれだけ大きいかを表す指標。
$\nu_1$ は、群間変動の自由度で、処理数（群の数）から 1 を引いた値。
$\nu_2$ は、群内変動の自由度で、データ総数から処理数を引いた値。
$n$ は、各群のデータ数（反復数）。
$x$ は、各データの測定値。
$\bar{x_i}$ は、i 番目の群の平均値。
$\bar{\bar{x}}$ は、すべてのデータの平均値（総平均）。

計算の流れ

群間変動を求める。各グループの平均と全体の平均の差の二乗和を計算する。

$群間変動 = 反復数 \times \sum(群平均 - 総平均) = n\sum(\bar{x_i}-\bar{\bar{x}})²$

群間変動の自由度を求める。処理数（グループ数）から 1 を引く。

$群間変動の自由度 = 処理数 - 1$

群内変動を求める。各データと、そのグループの平均の差の二乗和を計算する。

$群内変動\sum\sum(測定値 - 群平均)²$

群内変動の自由度を求める。総データ数から処理数を引く。

$群間変動の自由度 = データ総数 - 処理数$

要因分散を求める。群間変動を群間変動の自由度で割る。

$要因分散 = \frac{群間変動}{群間変動の自由度}$

誤差分散を求める。群内変動を群内変動の自由度で割る。

$誤差分散 = \frac{群内変動}{群内変動の自由度}$

F 値を求める。要因分散を誤差分散で割る。

$F値 = \frac{要因分散}{誤差分散}$

p 値を求めて、判定する。F 値と自由度を用いて、p 値を計算する。

R で実装する

対応のない場合を実装してみる。

デモデータの作成

3 つの水準（A, B, C）を持つ要因（group）と、2 つの水準（a, b）を持つ個体差（individual）によって分類された測定値（value）からなるデータを作る。

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

knitr::kable(t(sapply(levels(data$group), function(i){
  paste(subset(data, subset=group==i)$value,collapse=",")
})))

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

この表は、各群と各個体に対応する測定値を持つが、対応がない場合は、個体差(a, b)がいらない。このように、各水準が複数の値を持ち、水準以外の区別がないデータを「対応のある」と言う。

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

# 総平均
overall_mean <- mean(data$value)
# 群平均
group_means <- tapply(data$value, data$group, mean)

# 群間変動、「値-群平均」の偏差平方和
ss_between <- sum((group_means[data$group] - overall_mean)^2)
# 群間変動の自由度、水準数(処理数) - 1
df_between <- length(levels(data$group)) - 1
# 要因分散
ms_between <- ss_between / df_between

# 群内変動、誤差、「群平均-総平均」の偏差平方和
ss_within <- sum((data$value - group_means[data$group])^2)
# 群内変動の自由度、総データ数-水準数(処理数)
df_within <- length(data$value) - length(levels(data$group))
# 誤差分散
ms_within <- ss_within / df_within

# F値、要因分散/誤差分散
f_value <- ms_between / ms_within

# F値からp値を求める
p_value <- pf(f_value, df_between, df_within, lower.tail = FALSE)

# 結果の表示
cat("F-value:", round(f_value, 2), "\n")
cat("p-value:", round(p_value, 4), "\n")

F-value: 4.47
p-value: 0.0449

分散分析の関数を使う

# 対応のない一元配置分散分析
res <- aov(value ~ group, data = data)
summary(res)

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  13287    6644    4.47 0.0449 *
Residuals    9  13378    1486
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F 値は 4.47 で、p 値は 0.0449。

有意水準が 5%なら0.05 より小さいので、有意差あり、平均値に差がある、となる。

日々のつれづれ

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

分散分析のメモ・その１