はじめに

前回の記事で ¹ seabornにて公開されているpenguinsのデータセットを使って重回帰分析を行った。seabornにはpenguinsのデータセット以外にもいくつかのデータセットが公開されている。そこで今回はdiamondsのデータセットを使って一元配置分散分析を実施してみる。

その前に

まずは定義を確認しておきます。
一元配置分散分析：1つのカテゴリ変数の3つ以上のグループに基づいて、1つの連続従属変数の平均を比較。

必要なライブラリ、データのインポート

必要なパッケージを読み込み、データの確認を行います。

# パッケージのインポート
import pandas as pd
import seaborn as sns

# ダイヤモンドデータセットをseabonrから読み込む
diamonds = sns.load_dataset("diamonds", cache=False)

# 先頭10行を確認しておく。
print(diamonds.head(10))

出力例

   carat        cut color clarity  depth  table  price     x     y     z
0   0.23      Ideal     E     SI2   61.5   55.0    326  3.95  3.98  2.43
1   0.21    Premium     E     SI1   59.8   61.0    326  3.89  3.84  2.31
2   0.23       Good     E     VS1   56.9   65.0    327  4.05  4.07  2.31
3   0.29    Premium     I     VS2   62.4   58.0    334  4.20  4.23  2.63
4   0.31       Good     J     SI2   63.3   58.0    335  4.34  4.35  2.75
5   0.24  Very Good     J    VVS2   62.8   57.0    336  3.94  3.96  2.48
6   0.24  Very Good     I    VVS1   62.3   57.0    336  3.95  3.98  2.47
7   0.26  Very Good     H     SI1   61.9   55.0    337  4.07  4.11  2.53
8   0.22       Fair     E     VS2   65.1   61.0    337  3.87  3.78  2.49
9   0.23  Very Good     H     VS1   59.4   61.0    338  4.00  4.05  2.39

今回は各color（色）のグレードに対して、price（価格）の平均が統計的に有意な差があるかどうかを一元配置分散分析を使って確認していきたいと思います。

データの準備

color（色）とprice（価格）のデータを確認しておきます。
まずはcolor（色）から。

# 各色のグレードのダイヤモンドの数を確認
print(diamonds["color"].value_counts())

出力例

G    11292
E     9797
F     9542
H     8304
D     6775
I     5422
J     2808
Name: color, dtype: int64

次に、価格のヒストグラムを作成して分布を確認しておく。

import matplotlib.pyplot as plt
fig = sns.histplot(diamonds["price"])
fig.set_xlabel("Residual Value")
fig.set_title("Histogram of Residuals")
plt.show()

分布が大きく偏っていて高い値の長い裾が確認出来ます。そこで対数変換を行い歪んだ分布が対称的になるかを確認します。これは対数変換が正規性の仮定が必要な場合の統計分析とモデリングで特に役立つからですね。

# 価格の対数を取り、3列目に挿入します。
diamonds.insert(10, "log_price", [math.log(price) for price in diamonds["price"]])

# ヒストグラムを作成（省略）

正規分布に近づいてきました。よってprice（価格）ではなく、log_price（価格を対数化した値）を連続変数として扱うことにします。

探索的データ分析（EDA）

次に箱ひげ図を作成して基本的な探索的データ分析（EDA）を行います。

# 色のグレード別に価格の分布を示すボックス プロットを作成する
sns.boxplot(x = "color", y = "log_price", data = diamonds, hue="color", palette="Set1")
plt.show()

箱ひげ図より、ダイヤモンドの各色の log_price（価格を対数化した値）の分布に多くの重複していることが分かりました。しかし、統計的に有意に異なるかどうかはまだ分かりません。

分散分析（ANOVA：Analysis of Variance）を実行するには、回帰モデルを作成する必要があります。これを行う為に、statsmodels.api パッケージと ols() 関数を使用します。そして、変数Xがcolor（色）である、単純な線形回帰モデルを作成し、C() を使用してカテゴリとしてコード化します。次に、モデルをデータに適合させ、モデルの回帰分析の結果を取得します。

# statsmodels.api パッケージと ols() 関数をインポート
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 単純な線形回帰モデルを構築し、モデルを適合する
model = ols(formula = "log_price ~ C(color)", data = diamonds).fit()

# 分析結果の概要を取得
print(model.summary())

出力例

                            OLS Regression Results                            
==============================================================================
Dep. Variable:              log_price   R-squared:                       0.026
Model:                            OLS   Adj. R-squared:                  0.026
Method:                 Least Squares   F-statistic:                     237.8
Date:                Thu, 09 Jan 2025   Prob (F-statistic):          3.77e-301
Time:                        20:22:14   Log-Likelihood:                -76617.
No. Observations:               53940   AIC:                         1.532e+05
Df Residuals:                   53933   BIC:                         1.533e+05
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
=================================================================================
                    coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------
Intercept         7.6169      0.012    625.984      0.000       7.593       7.641
C(color)[T.E]    -0.0375      0.016     -2.370      0.018      -0.069      -0.006
C(color)[T.F]     0.1455      0.016      9.146      0.000       0.114       0.177
C(color)[T.G]     0.1727      0.015     11.219      0.000       0.143       0.203
C(color)[T.H]     0.3015      0.016     18.390      0.000       0.269       0.334
C(color)[T.I]     0.4061      0.018     22.250      0.000       0.370       0.442
C(color)[T.J]     0.5291      0.022     23.537      0.000       0.485       0.573
==============================================================================
Omnibus:                    11794.122   Durbin-Watson:                   0.059
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             2240.596
Skew:                           0.064   Prob(JB):                         0.00
Kurtosis:                       2.010   Cond. No.                         8.56
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

モデルの出力結果の P>|t| 列を確認すると、color（色）に関連付けられた係数βの p 値はすべて 0.05 未満であることが確認できます。つまり、color（色）は log_price（価格を対数化した値）に影響を及ぼす、と言えます。ただし、各color（色）間で価格の平均に大きな差（優位な差）がある、とは断定できません。そこで、一元配置分散分析を用います。

一元配置分散分析（One-way ANOVA）の実施

帰無仮説と対立仮説

まず、帰無仮説($H_0$)と対立仮説($H_1$)を確認しておきます。

帰無仮説：

H_0 :
𝑝𝑟𝑖𝑐𝑒_D = 𝑝𝑟𝑖𝑐𝑒_E = 𝑝𝑟𝑖𝑐𝑒_F = 𝑝𝑟𝑖𝑐𝑒_G = 𝑝𝑟𝑖𝑐𝑒_H = 𝑝𝑟𝑖𝑐𝑒_I = 𝑝𝑟𝑖𝑐𝑒_J

ダイヤモンドの価格の平均は色のグレードによって差はない。

対立仮説：

H_1 : Not (𝑝𝑟𝑖𝑐𝑒_D = 𝑝𝑟𝑖𝑐𝑒_E = 𝑝𝑟𝑖𝑐𝑒_F = 𝑝𝑟𝑖𝑐𝑒_G = 𝑝𝑟𝑖𝑐𝑒_H = 𝑝𝑟𝑖𝑐𝑒_I = 𝑝𝑟𝑖𝑐𝑒_J)

ダイヤモンドの価格の平均は色のグレードによって異なる。

Pythonで一元配置分散分析（One-way ANOVA）の実施

# 一元配置分散分析（One-way ANOVA）
print(sm.stats.anova_lm(model, typ = 1))

出力例

               df        sum_sq    mean_sq           F         PR(>F)
C(color)      6.0   1431.255783  238.54263  237.807767  3.767555e-301
Residual  53933.0  54099.661516    1.00309         NaN            NaN

PR(>F) 列の p 値が非常に小さいため、帰無仮説を棄却し、対立仮説を採用します。つまり、「ダイヤモンドの価格の平均は色のグレードによって異なる」が結論となります。

技術的な補足。

anova_lm の typ についてはここ ² とかここ ³ などを参照されたい。

Pythonを使って重回帰分析を行う https://qiita.com/tatsu_sekine/items/4c50a0b0e0257761cc2f ↩
Anova – Type I/II/III SS explained | R-bloggers https://www.r-bloggers.com/2011/03/anova-%e2%80%93-type-iiiiii-ss-explained/ ↩
r - How to interpret type I, type II, and type III ANOVA and MANOVA? - Cross Validated https://stats.stackexchange.com/questions/20452/how-to-interpret-type-i-type-ii-and-type-iii-anova-and-manova ↩