Pythonを用いた分散分析（ANOVA）あれこれ（1）

Pythonを用いて分散分析（ANOVA）を行う方法について勉強する機会がありましたので、分かったことをあれこれ纏めました。（1）では対応のない一元配置分散分析＋αについて纏めます。個人的な備忘録の色合いが強く、色々と端折っている部分があります。

※環境はJupyter Notebookの想定です。
※有意水準は全て5%（0.05）とします。

対応のない一元配置分散分析
多重比較法（Tukey-Kramer法）

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

使用するデータ

例えば100 gの樹脂サンプルを、3種の異なる溶媒A、B、Cにそれぞれ一定時間浸漬したときに、溶媒によって樹脂の溶け具合（浸漬後の重量）に差異が生じるか検証する実験を考えます。溶媒1種類につき5サンプルを評価し、浸漬後の重量を測定したところ、下表の結果が得られたとします。

	溶媒A	溶媒B	溶媒C
1	76	79	81
2	81	71	87
3	83	80	82
4	76	77	87
5	82	71	88

平均	79.6	75.6	84.8

表を見ると差異が生じていそうですが、分散分析できちんと検定してみます。この場合…

帰無仮説：全ての溶媒間において、浸漬後重量の平均値に差異は存在しない
対立仮説：少なくとも一組の溶媒間において、浸漬後重量の平均値に差異が存在する

となります。
　

Scipyのstatsモジュールを用いる方法

おそらく一番手軽な方法です。scipy.stats.f_oneway()関数を使用します。

import numpy as np
import pandas as pd
from scipy import stats

solventA = np.array([76, 81, 83, 76, 82])
solventB = np.array([79, 71, 80, 77, 71])
solventC = np.array([81, 87, 82, 86, 88])

#後々のために、一度データフレームにまとめておく
data = pd.DataFrame({
    'A' : solventA,
    'B' : solventB,
    'C' : solventC
})

fvalue, pvalue = stats.f_oneway(data['A'], data['B'], data['C'])
# *以下のようにNumPy配列を直接与えることも可能
#  fvalue, pvalue = stats.f_oneway(solventA, solventB, solventC)

print(f"F-val:{fvalue:.4},  P-val:{pvalue:.4}")

結果は以下の通りとなります。f_oneway関数が返すのはF値とp値のみであり、分散分析表のような形で群間変動や誤差変動などを得ることは出来ないようです。

F-val:8.0201,　P-val:0.006143

p値が5%（0.05）を下回ったので帰無仮説は棄却され、少なくとも一組の溶媒間において、浸漬後の重量の平均値には有意な差異が存在すると言えます。ただし具体的にどの溶媒間に差異が存在するのかは、分散分析では検定出来ません。
　

statsmodelsを用いる方法

統計解析用のライブラリであるstatsmodelsを用いて分散分析モデルを構築することが出来ます。まずはモデル化に向けて、先程のデータフレームを縦持ちに変換しておきます。

#先程のデータフレームをそのまま利用
data_melt = data.melt(var_name='Solvent', value_name='Weight')
data_melt

そしてstatsmodelsで、最小二乗法（ols）を用いて分散分析モデルを構築します。

from statsmodels.formula.api import ols
model = ols('Weight ~ Solvent', data=data_melt).fit()                
model.summary()

サマリーは以下のように表示されます。

後は構築したモデルを用いて、statsmodels.api.stats.anova_lm()により分散分析を実行出来ます。

import statsmodels.api as sm

#引数typは、平方和のタイプ指定（いわゆるANOVA Type I, II, IIIなどと分類されるもの）
#今回はtyp=2としたが、今回のように全ての溶媒に対して同数（5サンプル）のデータがある場合は
#結果に差異は生じない
anova = sm.stats.anova_lm(model, typ=2)
anova

結果として、分散分析表がpandasデータフレームの形で得られます。F値、p値はScipyを用いたときと一致しました。また群間変動や誤差変動（sum_sq）、自由度（df）も出力されています。

　

補足：分散分析モデルの構築について

最初は、何故statsmodelsでは最小二乗法（ols）を用いたモデル化が必要になるのか？について理解が出来ていませんでしたが、色々と調べるうちに背景が分かってきました。数学的な正確性や厳密性には自信がありませんが、自分なりに理解出来たことをメモしてみます（以下の外部サイトが参考になりました）
一般化線形モデル（PDF資料）
正規方程式の導出と計算例 | 高校数学の美しい物語
 pythonで統計学基礎：03 検定・分散分析 | 化学の新しいカタチ

今回の場合、分散分析モデルは以下のように定式化出来ます。溶媒Aにおける浸漬後の重量データを基準として、溶媒B、溶媒Cの場合はそれぞれ補正が加わるようなイメージです。

- $y$ ：浸漬後の重量
- $x_1$ ：溶媒Bの場合は1、そうでない場合は0となる変数
- $x_2$ ：溶媒Cの場合は1、そうでない場合は0となる変数
- $e$ ：誤差項
- $\beta_0, \beta_1, \beta_2$ ：パラメータ

ここで溶媒Xを用いたときのサンプルiについて、浸漬後の重量を $y_{Xi}$ 、誤差項を $e_{Xi}$ と表現すると、以下のように行列で表記が出来ます。

$\begin{bmatrix} y_{A1} \\ y_{A2} \\ y_{A3} \\ y_{A4} \\ y_{A5} \\ y_{B1} \\ y_{B2} \\ y_{B3} \\ y_{B4} \\ y_{B5} \\ y_{C1} \\ y_{C2} \\ y_{C3} \\ y_{C4} \\ y_{C5} \end{bmatrix}　=　\begin{bmatrix} 1 & 0 & 0\\ 1 & 0 & 0\\ 1 & 0 & 0\\ 1 & 0 & 0\\ 1 & 0 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & 0 & 1\\ \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix} + \begin{bmatrix} e_{A1} \\ e_{A2} \\ e_{A3} \\ e_{A4} \\ e_{A5} \\ e_{B1} \\ e_{B2} \\ e_{B3} \\ e_{B4} \\ e_{B5} \\ e_{C1} \\ e_{C2} \\ e_{C3} \\ e_{C4} \\ e_{C5} \end{bmatrix}$

これをまとめて、以下のように表記することが出来ます。

$\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E}$

これは一般線形モデルであり、誤差の二乗和が最小となる時のパラメータ $\mathbf{B}$ （ $\mathbf{B} = [\beta_0 \ \beta_1 \ \beta_2]^{T}$ ）は、以下のように正規方程式を解くことで求めることが出来ます。

$\mathbf{B} = \left( \mathbf{X}^{T} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \mathbf{Y}$

条件として $\left( \mathbf{X}^{T} \mathbf{X} \right)$ が正則である（逆行列を持つ）必要がありますが、今回のような定式化であれば満たすことが出来ます。実際にこの計算を自力でやってみました。

#先程の縦持ち化したデータフレームから、浸漬後の重量の値を取得
Y = data_melt['Weight'].to_numpy()

X = np.array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], 
                [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0], 
                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1]]).T

#正規方程式を解く
B = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print(B)

結果は以下の通りであり、 $\beta_0=79.6, \ \beta_1=-4.0, \ \beta_2=5.2$ となります。

[79.6 -4. 5.2]

これは先程、statsmodelsのolsが出力したサマリーにおける係数（coef）の値とピタリ一致します。つまり、olsを用いた分散分析モデル構築の中身は、このような演算ということになります。

$\beta_0=79.6$ は溶媒Aにおける浸漬後の重量の平均そのものであり、そこに $\beta_1=-4.0$ を加えた値は溶媒Bにおける平均（75.6）、 $\beta_1=5.2$ を加えた値は溶媒Cにおける平均（84.8）となります。
　

多重比較法（Tukey-Kramer法）

先述の通り分散分析では、具体的にどの溶媒間に差異が存在するのかは検定することが出来ません。そのような目的のためには多重比較法を用いる必要があります。

ただし「事前に分散分析を行い、その結果が有意であった場合に事後検定として多重比較法を行い、具体的にどの群間に有意差があるのかを検定する」という考え方は正しくないようです。

多重比較法には様々な種類があり、確かに分散分析を事前に行うことを前提とした手法（Scheffe法など）もありますが、一方で分散分析とは独立して群間の有意差の検定に適用可能な手法（Bonferroni法、Tukey-Kramer法など）もあります。もちろん後者においては、分散分析と結果が整合しない場合（分散分析では有意なのに、多重比較法では群間の有意差が出てこない、あるいはその逆）もあり得るようです。

分散分析と多重比較法の関係については、以下の外部サイトが参考になりました。
分散分析の下位に多重検定を置くな
 検定の多重性 | ブログ | 統計WEB
統計検定を理解せずに使っている人のために III（PDF資料）

ここでは溶媒への樹脂浸漬データに対して、Tukey-Kramer法を用いた多重比較を行なってみます。

#先程の縦持ち化したデータフレームを使用
from statsmodels.stats.multicomp import pairwise_tukeyhsd
pairwise_tukeyhsd(data_melt["Weight"], data_melt["Solvent"]).summary()

サマリーは以下のように表示されます。溶媒B、C間では、p値（p-adj）が5%（0.05）を下回っており、浸漬後重量の平均値に有意差ありとなります（その場合「reject」がTrueとなります）。

今回の結果は分散分析とも整合しておりますが、先述の通りTukey-Kramer法は分散分析とは独立した手法であり、必ず整合するとは限らないことには注意が必要です。

工場裏のアーカイブス

素人によるiPhoneアプリ開発の学習記　あと機械学習とかM5Stackとか