勉強のためのメモ。2025年6月25日更新。
Contents
データ
#R内での上のデータの生成
N=15
set.seed(-6152847)
a <- rnorm(N, mean=5.080351, sd=0.1375497)
set.seed(3448258)
b <- rnorm(N, mean=5.185921, sd=0.05358547)
set.seed(-9846648)
c <- rnorm(N, mean=5.186056, sd=0.05089308)
set.seed(-1377657)
d <- rnorm(N, mean=5.186409, sd=0.06550713)
set.seed(7612155)
e <- rnorm(N, mean=5.238295, sd=0.07842281)
set.seed(NULL)
value <- c(a, b, c, d, e)
value <- round(value*10, digit=1)
group <- factor(c(rep("A", N), rep("B", N), rep("C", N), rep("D", N), rep("E", N)))
| 群 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 19 | 11 | 12 | 13 | 14 | 15 |
| A | 51.3 | 50.7 | 53.7 | 50.6 | 49.0 | 50.7 | 53.7 | 49.8 | 51.1 | 51.3 | 52.2 | 50.8 | 49.0 | 51.8 | 51.0 |
| B | 52.6 | 51.0 | 52.5 | 51.4 | 52.3 | 52.5 | 52.0 | 51.2 | 52.2 | 51.5 | 51.5 | 51.9 | 52.3 | 51.4 | 51.1 |
| C | 52.2 | 52.4 | 52.4 | 51.2 | 51.6 | 52.6 | 51.6 | 52.2 | 51.1 | 51.8 | 52.8 | 52.1 | 51.6 | 51.6 | 52.3 |
| D | 53.3 | 51.4 | 51.8 | 52.7 | 51.3 | 51.8 | 51.1 | 52.2 | 50.8 | 51.2 | 52.8 | 52.1 | 52.7 | 53.0 | 51.6 |
| E | 53.9 | 52.5 | 52.7 | 52.1 | 52.6 | 54.4 | 52.3 | 52.0 | 52.5 | 52.6 | 51.3 | 51.4 | 52.7 | 51.5 | 51.7 |
| 群 | サイズ \(n_i\) | 計 \(T_i\) | 平均 \(\bar{x}_i\) | 不偏分散 \(V_i\) |
| A | 15 | 766.7 | 51.11333 | 1.891238 |
| B | 15 | 777.4 | 51.82667 | 0.3078095 |
| C | 15 | 779.5 | 51.96667 | 0.2580952 |
| D | 15 | 779.8 | 51.98667 | 0.5955238 |
| E | 15 | 786.2 | 52.41333 | 0.7355238 |
全データ数\(N = 75\)。データの総和\(T = 3889.6\)。全データの平均 \(\bar{\bar{x}} = 51.86133\)。群の数\(k = 5\)
等分散を仮定する分散分析の実行
分散分析の根幹は、総平方和が群間のばらつき（群間平方和）と標本誤差のばらつき（群内平方和）に分解できるという事実。
用語
- 総平方和（\(S_T\)） — それぞれの値と全体の平均値との差を二乗して足し合わせたもの。
- 群間平方和（\(S_A\)、A間平方和） — 群内平均と全体平均の差を二乗して足し合わせたもの。
- 群内平方和（\(S_E\)、誤差平方和） — それぞれの値と各群内の平均値との差を二乗して足し合わせたもの。
各平方和の計算
$$ CT = \frac{T^2}{N} = 201719.8$$
$$S_T = \sum_{i=1}^k{} \sum_{j=1}^{n_i}(x_{ij}-\bar{\bar{x}})^2 = \sum_{i=1}^k{} \sum_{j=1}^{n_i} x^2_{ij} - CT = 66.41787$$
$$S_A = \sum_{i=1}^k{} \sum_{j=1}^{n_i}(\bar{x}_{i}-\bar{\bar{x}})^2 = \sum_{i=1}^k{}\frac{T_i^2}{n_i} - CT = 13.3832 $$
$$S_B = \sum_{i=1}^k{} \sum_{j=1}^{n_i}(x_{ij}-\bar{x})^2 = S_T - S_A = 53.03467$$
分散分析表
| 要因 | 平方和 \(S\) | 自由度 \(\phi\) | 平均平方 \(V\) | \(F_0\) |
| \(A\) | 13.3832 | 4 | \(S_A/\phi_A =\) 3.3458 | \(V_A/V_E\) = 4.416093 |
| \(E\) | 53.03467 | 70 | \(S_E/\phi_E =\) 0.7576381 | |
| 計 | 66.41787 | 74 |
データを平均0、分散1に標準化 (standardization) するには、値と平均の差を取って (\(X - \bar{x}\))、標準偏差 (\(s\)) で割ります。このスケーリングされた値\(\frac{X-\bar{x}}{s}\) が正規分布（平均0、分散1）からサンプリングされているとすると、その平方和（\(\sum_{i=1}^k \frac{(x_i - \bar{x})^2}{s^2}\)）は自由度 \(k-1\) の\(\chi^2\)分布（カイじじょうぶんぷ）に従います。正規分布からサンプリングされた値の平方和をそれらの自由度\(k_1, k_2\)で割った値の比は自由度 (\(k_1, k_2\)) のF分布に従います。
したがって、平方和\(S_A, S_E\)（を分散\(s^2\)で割った値）はそれぞれ自由度\(\phi_A, \phi_E\) の\(\chi^2\)分布に従い、F値は各群の平均値の不偏分散（\(S_A/\phi_A\)）と誤差の不偏分散（\(S_E/\phi_E\)）の比$$F = \frac{\frac{S_A}{\phi_A}}{\frac{S_E}{\phi_E}}$$として計算できます。分散\(s^2\)は比を取った時に打ち消されるので計算式には現われていません。
Rでp値を計算します。
> pf(4.416093, 4, 70, lower.tail=FALSE)
[1] 0.003058812
> oneway.test(value ~ group, var.equal=TRUE) # デフォルトはvar.equal=FALSEで等分散を仮定しない
One-way analysis of means
data: value_ and group
F = 4.4161, num df = 4, denom df = 70, p-value = 0.003059
> anova(aov(value ~ group)) # aov関数は等分散を仮定する分散分析しか実行できない
Analysis of Variance Table
Response: value
Df Sum Sq Mean Sq F value Pr(>F)
group 4 13.383 3.3458 4.4161 0.003059 **
Residuals 70 53.035 0.7576
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
等分散を仮定しない分散分析の実行
$$w_i = \frac{n_i}{V_i}$$
$$\hat{x} = \frac{\sum_i^k w_i \bar{x}_i}{\sum_i^k w_i}$$
とします。母平均が全て等しく、標本サイズ（\(n\)）が大きい（すなわち群内自由度が大きい）時、統計量$$\sum_i^k w_i(\bar{x}_i - \hat{x})^2$$は自由度\(k-1\)の\(\chi^2
\)分布に従います。
\(k=2\)の時
\(k=2\)の時は
$$\sum_{i=1}^k w_i(\bar{x}_i - \hat{x})^2 = v^2 = \left( \frac{\bar{x}_1 - \bar{x}_2}{\sqrt{\frac{V_1}{n_1}+\frac{V_2}{n_2}}} \right)^2$$
2つの群の母平均が等しい時は、\(v\)の分子は標準偏差\(\sqrt{\sigma_1^2/n_1 + \sigma_2^2/n_2} \)を持つ正規偏差です。\(v\)の分母は\(\sigma^2\)を不偏分散\(s^2\)で置き変えることによって得られるこの標準偏差の推定値です。
\(k>2\)の時
\(k>2\)の時は統計量の分布を近似しなければいけません。
$$v^2 = \frac{\sum_{i=1}^k \frac{w_i(\bar{x}_i-\hat{x})^2}{k-1}}{\left[ 1 + \frac{2(k-2)}{k^2-1}\sum_{i=1}^k\frac{1}{n_i-1}\left(1-\frac{w_i}{\sum_{i=1}^k w_i}\right)^2 \right]}$$
自由度は
$$\nu_1 = (k-1)$$
$$\nu_2 = \left[ \frac{3}{(k^2-1)}\sum_{i=1}^k\frac{1}{n_i-1}\left(1- \frac{w_i}{\sum_{i=1}^k w_i}\right)^2 \right]^{-1}$$
計算
| 群 | サイズ \(n_i\) | 平均 \(\bar{x}_i\) | 不偏分散 \(V_i\) | 平均値の推定分散 \(V_i/n_i\) |
\(w_i\) (\(n_i/V_i)\) |
| A | 15 | 51.1333 | 1.891238 | 0.12608254 | 7.931312 |
| B | 15 | 51.82667 | 0.3078095 | 0.02052063 | 48.731436 |
| C | 15 | 51.96667 | 0.2580952 | 0.01720635 | 58.118081 |
| D | 15 | 51.98667 | 0.5955238 | 0.0390159 | 25.187910 |
| E | 15 | 52.41333 | 0.7355238 | 0.04903492 | 20.393629 |
| 計 | 160.3624 |
もし平均値がすべて等しいならば、全体平均の最良の推定値は\(\hat{x} = (\sum_{i=1}^k w_i\bar{x}_i)(\sum_{i=1}^k w_i) = 51.94186\)となります。統計量の手計算に必要な値を下の表に示します。任意の原点として\(x =52\)が使用されています。
| 群 | \(w_i(\bar{x}_i-52)\) | \(w_i(\bar{x}_i-52)^2\) | \(w_i/\sum w_i\) | \((1- w_i/\sum w_i)^2\) |
\(\frac{(1- w_i/\sum w_i)^2}{n_i-1}\) |
| A | -7.0324303 | 6.235421493 | 0.04945869 | 0.9035288 | 0.06453777 |
| B | -8.4467822 | 1.464108911 | 0.30388324 | 0.4845785 | 0.03461275 |
| C | -1.9372694 | 0.064575646 | 0.36241720 | 0.4065118 | 0.02903656 |
| D | -0.3358388 | 0.004477851 | 0.15706871 | 0.7105332 | 0.05075237 |
| E | 8.4293668 | 3.484138288 | 0.12717216 | 0.7618284 | 0.05441632 |
| 計 | -9.322954 | 11.25272 | 0.2333558 |
$$\sum_{i=1}^k w_i(\bar{x}_i - \hat{x})^2 = 11.25272 - (9.322954)^2/160.3624 = 10.71071$$
\(v^2\)を求める式に代入すると、
$$ v^2 = \frac{1/(k - 1) (10.71071) }{1 + \frac{2(k-2)}{k^2-1}(0.2333558) } = \frac{\frac{1}{4}(10.71071)}{1+ \frac{6}{24}(0.2333558)}=2.530076$$
$$\nu_2 = \left[\frac{3}{24}(0.2333558)\right]^{-1} = \frac{8}{0.2333558} = 34.28241$$
Rでp値を計算します。
> pf(2.530076, 4, 34.28241, lower.tail=FALSE) [1] 0.05823581 > oneway.test(value ~ group) One-way analysis of means (not assuming equal variances) data: value_ and group F = 2.5301, num df = 4.000, denom df = 34.282, p-value = 0.05824
（了）