分散分析の概要

以下のような問題を考えます。

【例題】

因子Aについて4つの水準A1~A4があり、それぞれの水準で以下のような観測が得られた。
水準間に違いはあると言えるか、有意水準5%で検定せよ。

level \ n	1	2	3	4	5
A1	15	13	15	16	14
A2	18	17	16	15	18
A3	19	16	17	18
A4	17	15	16

2標本の平均値の差の検定にはt検定が用いられますが、上記の問題のように3群以上からなるデータの場合には「分散分析」を行います。この場合、因子がAのみであるため「一元配置分散分析」です。

帰無仮説$H_0$は「各水準の母平均は等しい」、対立仮説$H_0$「少なくとも一つの水準について母平均が異なる」となります。

分散分析の実装

パッケージを利用して分散分析を行う場合には statsmodels.stats.anova などが用いられますが、今回はパッケージを利用せずに実施します。全体の流れとしては次の通りです。

Step1. 総平方和を求める
Step2. 水準間平方和、誤差平方和を求める
Step3. F統計量を求める
Step4. F分布から棄却域を決定する
Step5. 帰無仮説を棄却するかどうか判断する

import numpy as np
from scipy.stats import f

Y = [
    [15, 13, 15, 16, 14], 
    [18, 17, 16, 15, 18], 
    [19, 16, 17, 18], 
    [17, 15, 16]
]

Step1. 総平方和を求める

まず、データ全体のデータの個数と平均値を求めます。

n_instances = 0
total = 0
for y_A in Y:
    n_instances += len(y_A)
    total += sum(y_A)

total_mean = total / n_instances
print('データの個数: {}\nデータ全体の平均値: {}'.format(n_instances, total_mean))

Out

データの個数: 17
データ全体の平均値: 16.176470588235293

応答yについて、全体の平均値からの偏差について二乗和をとり、総平方和$S_T$を次のように求めます。

$$S_T=\sum_{i=1}^{a}\sum_{j=1}^{n}(y_{ij}-\bar{y})^2$$

S_T = 0
for i, y_A in enumerate(Y):
    for j in range(len(y_A)):
        S_T += (Y[i][j] - total_mean)**2

print('S_T = {}'.format(S_T))

Out

S_T = 40.47058823529411

Step2. 水準間平方和、誤差平方和を求める

続いて水準間平方和$S_A$を次のように求めます。水準ごとの平均のばらつきが大きいほど$S_A$が大きくなります。

$$S_A=\sum_{i=1}^{a}\sum_{j=1}^{n}(\bar{y}_{A_i}-\bar{y})^2$$

means = []
for y_A in Y:
    mean_A = sum(y_A)/len(y_A)
    means.append(mean_A)

S_A = 0
for i, mean_A in enumerate(means):
    S_A += (mean_A - total_mean)**2 * len(data[i])

print('S_A = {}'.format(S_A))

Out

S_A = 21.47058823529413

総平方和$S_T$から水準間平方和$S_A$を引くと、誤差平方和$S_E$に等しくなります。

S_E = S_T - S_A
print('S_E = {}'.format(S_E))

Out

S_E = 18.99999999999998

Step3. F統計量を求める

総平方和, 水準間平方和, 誤差平方和の自由度$\phi_T, \phi_A, \phi_E$はそれぞれ以下のように求められます。

phi_T = n_instances - 1
phi_A = len(Y) - 1
phi_E = phi_T - phi_A

print(phi_T, phi_A, phi_E)

Out

16 3 13

平方和$S_A, S_E$を自由度$\phi_A, \phi_E$で割り、水準間の変動$V_A=S_A/\phi_A$、および誤差による変動$V_E=S_E/\phi_E$を求めます。

帰無仮説が成り立つ場合には、次のF統計量は自由度 ($\phi_A, \phi_E$) のF分布に従います。

$$F=\frac{V_A}{V_E}=\frac{S_A/\phi_A}{S_E/\phi_E}$$

F = V_A / V_E
print('F = {}'.format(F))

Out

F = 4.896800825593403

Step4. F分布から棄却域を決定する

自由度 (3, 13) のF分布の上側5%点を求めます。

f_distribution = f(phi_A, phi_E)
f_distribution.ppf(0.95)

Out

3.4105336446278476

よって棄却域は $F \geq 3.411$ となります。

Step5. 帰無仮説を棄却するかどうか判断する

Step3で求めたF統計量は棄却域に含まれるので、帰無仮説は棄却されます。
すなわち、「水準間の母平均には違いがあるといえる」と結論づけられます。

Pythonで一元配置分散分析