分散分析（ANOVA）は、3つ以上のグループの平均値に統計的に有意な差があるかどうかを検定する手法です。通常のt検定を繰り返し行うと、第一種の過誤（偽陽性）の確率が累積して増大する「多重比較の問題」が発生します。分散分析を用いることで、このリスクを抑えつつ、一度の検定で全体としての有意差を判断できます。

前提条件と基礎概念

分散分析を適用する際は、データが正規分布に従うこと、および各グループの分散が等しいこと（等分散性）が前提となります。分析では以下の概念が重要です。

• F値：群間分散（要因の効果）を群内分散（誤差）で割った値です。F値が大きいほど、要因による影響が誤差に比べて大きいことを示します。
• p値：F分布の累積分布関数を用いて計算され、通常0.05未満であれば帰無仮説（すべての平均値が等しい）を棄却します。

実装環境の準備

数値計算、データ加工、可視化、統計モデリングに必要なライブラリを読み込みます。

import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf

sns.set_theme()

データの準備と可視化

天候がビールの売上に与える影響を検証するためのサンプルデータを作成します。サンプルサイズが小さいため、分布を確認しやすいよう箱ひげ図を使用します。

# データセットの作成
conditions = ['Cloudy', 'Cloudy', 'Rainy', 'Rainy', 'Sunny', 'Sunny']
sales = [6, 8, 2, 4, 10, 12]
df = pd.DataFrame({'Sales': sales, 'Weather': conditions})

# データの可視化
sns.boxplot(x='Weather', y='Sales', data=df)
plt.show()

# 各グループの平均値の確認
print(df.groupby('Weather').mean())

分散分析の手動計算

統計ライブラリを使わずに、F値とp値を導出するプロセスを確認します。これにより、分散分析の内部構造を理解できます。

1. 平方和の計算

各グループの平均値を予測値とみなし、全平均からの偏差平方和を「群間平方和」、実測値と予測値の差（残差）の平方和を「群内平方和」として計算します。

# 全体の平均値
grand_mean = df['Sales'].mean()

# 各天候ごとの平均値（予測値）
group_means = df.groupby('Weather')['Sales'].transform('mean')

# 群間平方和
ss_between = np.sum((group_means - grand_mean) ** 2)

# 群内平方和
ss_within = np.sum((df['Sales'] - group_means) ** 2)

print(f"群間平方和: {ss_between}")
print(f"群内平方和: {ss_within}")

2. 自由度と分散の計算

群間および群内の自由度を求め、それぞれの平均平方（分散）を計算します。

# パラメータ設定
k = df['Weather'].nunique()  # グループ数
N = len(df)                  # 全データ数

# 自由度
df_between = k - 1
df_within = N - k

# 平均平方（分散）
ms_between = ss_between / df_between
ms_within = ss_within / df_within

print(f"群間分散: {ms_between}")
print(f"群内分散: {ms_within}")

3. F値とp値の算出

F値は群間分散と群内分散の比率です。この値に基づき、F分布からp値を求めます。

f_value = ms_between / ms_within
p_value = 1 - stats.f.cdf(f_value, df_between, df_within)

print(f"F値: {f_value}")
print(f"p値: {p_value}")

計算結果としてp値が0.05を下回れば、天候が売上に有意な影響を与えていると判断できます。ただし、実務ではより多くのサンプルサイズが必要です。

線形モデルとしての解釈

分散分析は、説明変数がカテゴリカル変数である線形回帰モデルとして捉えることができます。例えば、天候が晴れかどうかを表すダミー変数（0または1）を用いることで、各天候の影響度を係数として推定します。

Statsmodelsを用いた実装

Pythonの統計ライブラリを使用すれば、これらの手順を自動化できます。

分散分析表の出力

`statsmodels`の`ols`関数でモデルを作成し、`anova_lm`で分散分析表を出力します。

# モデルの定義と適合
model = smf.ols('Sales ~ Weather', data=df).fit()

# 分散分析表の表示
anova_result = sm.stats.anova_lm(model, typ=2)
print(anova_result)

モデル係数の確認

モデルの係数を確認することで、各条件ごとの予測値がどのように構成されているか理解できます。

# 係数の出力
print(model.params)

# 予測値と残差
print("予測値:")
print(model.fittedvalues)
print("残差:")
print(model.resid)

このモデルでは、基準となるカテゴリ（通常はアルファベット順で最初のもの）の平均が切片となり、他のカテゴリの係数はその基準との差分として表されます。

回帰モデルへの応用

分散分析の考え方は、説明変数が連続変数である回帰モデルにも適用できます。この場合、群間変動は「回帰変動」、群内変動は「残差変動」と呼ばれます。

# 連続変数を用いたデータセット
df_temp = pd.DataFrame({
    'Sales': np.array([45.3, 59.3, 40.4, 38.0, 37.0, 40.9, 60.2, 63.3, 51.1, 44.9, 47.0,
                       53.2, 43.5, 53.2, 37.4, 59.9, 41.5, 75.1, 55.6, 57.2, 46.5, 35.8,
                       51.9, 38.2, 66.0, 55.3, 55.3, 43.3, 70.5, 38.8]),
    'Temperature': np.array([20.5, 25.0, 10.0, 26.9, 15.8, 4.2, 13.5, 26.0, 23.3, 8.5, 26.2,
                             19.1, 24.3, 23.3, 8.4, 23.5, 13.9, 35.5, 27.2, 20.5, 10.2, 20.5,
                             21.6, 7.9, 42.2, 23.9, 36.9, 8.9, 36.4, 6.4])
})

# 線形回帰モデルの適用
reg_model = smf.ols('Sales ~ Temperature', data=df_temp).fit()

# 分散分析表の出力
print(sm.stats.anova_lm(reg_model, typ=2))

この結果から、気温が売上に対して有意な線形関係を持っているかどうかを判断できます。