日々のつれづれ

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

分散分析のメモ・その2

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

一元配置分散分析

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

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

対応がある場合の一元配置分散分析

1 つの要因につき同じ個体や対象を複数回測定したとき「対応がある」と言う。

各群のデータが同じ個体や対象ごとに対応関係があるので、同じ人に異なる処置を施してその効果を測定した場合や、同じ物体に異なる条件で実験した場合などが該当する。

F 値の計算では、この個体や対象の違いによる変動(個体間変動)がばらつきに加わるので、群内変動からこれを除いた調整済み誤差分散を使う。

F 値の式

\begin{eqnarray}

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

\end{eqnarray}
  • 調整済み群内変動は、群内変動から個体間変動を引いた値。要因とは別の影響によるデータのばらつきを示す。
  • $\nu_1$ は、群間変動の自由度で、処理数(群の数)から 1 を引いた値。
  • $\nu_2$ は、調整済み群内変動の自由度で、群内変動の自由度から個体間変動の自由度を引いた値。
  • $n$ は、各群のデータ数(反復数)。
  • $\bar{x_i}$ は、i 番目の群の平均値。
  • $\bar{\bar{x}}$ は、すべてのデータの平均値(総平均)。
  • $x$ は、各データの測定値。
  • $\bar{x_j}$ は、j 番目の個体や対象の平均値。

計算の流れ

  1. 群間変動を求める。各グループの平均と全体の平均の差の二乗和を計算する。
  2. 群間変動の自由度を求める。処理数(グループ数)から 1 を引く。
  3. 群内変動を求める。各データと、そのグループの平均の差の二乗和を計算する。
  4. 群内変動の自由度を求める。総データ数から処理数を引く。
  5. 個体間変動を求める。各個体の平均と全体の平均の差の二乗和を計算する。 $個体間変動 = \sum\sum(個体別平均 - 総平均)2 = \sum\sum(\bar{x_j}-\bar{\bar{x}})2$
  6. 調整済み郡内変動を求める。郡内変動から個体間変動を引く。 $調整済み郡内変動 = 郡内変動 - 個体間変動$
  7. 個体間変動の自由度を求める。個体数から 1 を引く。 $個体間変動の自由度 = 個体数 - 1$
  8. 調整済み郡内変動の自由度を求める。郡内変動の自由度から個体間変動の自由度を引く。 $調整済み郡内変動の自由度 = 群内変動の自由度-個体間変動の自由度 = \nu_{2}$
  9. 調整済み誤差分散を求める。調整済み郡内変動を調整済み郡内変動の自由度で割る。 $調整済み誤差分散 = \frac{調整済み群内変動}{調整済み群内変動の自由度}$
  10. F 値を求める。要因分散を調整済み誤差分散で割る。 $F_{(\nu_1,\nu_2)} = \frac{要因分散}{調整済み誤差分散}$
  11. p 値を求めて、判定する。F 値と自由度を用いて、p 値を計算する。

R で実装する

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

デモデータの作成

デモ用のデータを作成します。このデータは、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(sapply(levels(data$group), function(i) {
  sapply(levels(data$individual), function(j) {
    paste(subset(data, subset = group == i & individual == j)$value, collapse = ",")
  })
}))
A B C
a 31,15 37,14 90,91
b 14,3 118,153 91,92

それぞれの個体 (a と b) が、それぞれ要因の各水準 (A、B、C) の値を持っている。このように、1 つの個体が複数の値を持つデータを「対応のある」と言う。

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

# 総平均
overall_mean <- mean(data$value)
# 個体間平均
individual_mean <- tapply(data$value, data$individual, mean)
# 個体間変動、「個体別平均-総平均」の偏差平方和
ss_individual = sum((individual_mean[data$individual] - overall_mean)^2)
# 個体間変動の自由度、水準数(個体数) - 1
df_individual <- length(levels(data$individual)) - 1

# 個体差調整済み郡内変動、郡内変動-個体間変動
ss_within_new <- ss_within - ss_individual
# 個体差調整済み郡内変動の自由度、郡内変動の自由度-個体間変動の自由度=(総データ数-処理数)-(個体数-1)
df_within_new <- df_within - df_individual
# 個体差調整済み誤差分散
ms_within_new <- ss_within_new / df_within_new

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

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

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

分散分析の関数を使う

対応のある場合は、要因と個体を$+$で繋ぐ。

# 対応のある一元配置分散分析
res <- aov(value ~ group + individual, data = data)
summary(res)
            Df Sum Sq Mean Sq F value Pr(>F)
group        2  13287    6644   5.173 0.0362 *
individual   1   3104    3104   2.417 0.1586
Residuals    8  10274    1284
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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