分散分析の基礎 | Logics of Blue

https://logics-of-blue.com/anova-foundation/

前準備

In [1]:
# ライブラリの読み込み
import numpy as np
import pandas as pd
import scipy as sp
import statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf
In [2]:
# ライブラリのバージョン
print('np ',np.__version__)
print('pd ',pd.__version__)
print('sp ',sp.__version__)
print('sm ',statsmodels.__version__)

np  1.18.5
pd  1.0.5
sp  1.5.0
sm  0.11.1

分散分析の実行

In [3]:
# データの読み込み
data = pd.read_csv('data.csv')
print(data)

   weight food
0       2    A
1       4    A
2      10    B
3      12    B
4       6    C
5       8    C
In [4]:
# モデル化
anova_mod = smf.ols('weight ~ food', data=data).fit()

# 分散分析表
anova_table = sm.stats.anova_lm(anova_mod, type=2)
print(anova_table)

           df  sum_sq  mean_sq     F    PR(>F)
food      2.0    64.0     32.0  16.0  0.025095
Residual  3.0     6.0      2.0   NaN       NaN
In [5]:
# F分布を使ってp値を計算する
1 - sp.stats.f.cdf(16, 2, 3)
Out[5]:
0.02509457330439091
In [6]:
# 分散分析表はDataFrameなので結果の取り出しは容易
type(anova_table)
Out[6]:
pandas.core.frame.DataFrame
In [7]:
# F比
anova_table['F']
Out[7]:
food        16.0
Residual     NaN
Name: F, dtype: float64
In [8]:
anova_table.loc['food', 'F']
Out[8]:
16.0
In [9]:
# モデルから直接取り出してもよい
anova_mod.fvalue
Out[9]:
16.0

PB検定

In [10]:
# 帰無仮説が正しい(foodによって体重が変わらない)ことを仮定したモデル
anova_null = smf.ols('weight ~ 1', data=data).fit()
In [11]:
# 総平均
anova_null.params
Out[11]:
Intercept    7.0
dtype: float64
In [12]:
np.mean(data.weight)
Out[12]:
7.0
In [13]:
# モデルの標準偏差
np.sqrt(anova_null.scale)
Out[13]:
3.7416573867739413
In [14]:
np.std(data.weight, ddof=1)
Out[14]:
3.7416573867739413
In [15]:
# 平均7、標準偏差3.741657の正規分布に従う確率変数を生成
np.random.seed(1)
mu = np.mean(data.weight)
sd = np.std(data.weight, ddof=1) 
size = len(data)

norm_instance = sp.stats.norm(loc=mu, scale=sd)
res_rnorm = norm_instance.rvs(size=size)

print(res_rnorm)

[13.07774383  4.7110171   5.02376226  2.98531903 10.23805885 -1.61156927]
In [16]:
# このシミュレーションデータに対してF比を計算すると、
# 小さなF比が得られやすい

# DataFrameに直す(statsmodels.regression.linear_model.OLSを使う方法もある)
sim_df = pd.DataFrame({
    'weight': res_rnorm,
    'food': data.food
})

# 分散分析の実行
mod_sim = smf.ols('weight ~ food', data=sim_df).fit()
print('F比', mod_sim.fvalue)
print(sm.stats.anova_lm(mod_sim, type=2))

F比 0.4193728369926777
           df      sum_sq    mean_sq         F    PR(>F)
food      2.0   29.995090  14.997545  0.419373  0.690872
Residual  3.0  107.285527  35.761842       NaN       NaN
In [17]:
# PB検定の実行
# 時間がかかるので注意

n_sim = 50000                 # シミュレーションする回数
f_ratio_vec = np.zeros(n_sim) # F比を保管する容れ物
exp_vec = data.food           # 水準

np.random.seed(1)
for i in range(0, n_sim):
    # シミュレーションにより体重データを生成
    # このデータは、体重が餌によって変化しないことを想定している
    simlated_weight = norm_instance.rvs(size=size)
    sim_df = pd.DataFrame({
        'weight': simlated_weight,
        'food': exp_vec
    })
    # モデル化と分散分析表の出力
    mod_sim = smf.ols('weight ~ food', data=sim_df).fit()
    # F比を保管
    f_ratio_vec[i] = mod_sim.fvalue
In [18]:
# F比が16を超えた割合
sum(f_ratio_vec > 16) / n_sim
Out[18]:
0.02458