Pythonでanova
この本をPythonで書き直すシリーズです。(クリックするとアマゾンに飛びます)
後半戦が重くてなかなか進んでません。
（そろそろ本格的に統計をやらなければいけない気がしてきた）
今回はanova(analysis of variance:分散分析)です。
あまりPythonでやる日本語情報は見つかりませんでした。
モジュール読み込み
import pandas as pd import statsmodels.api as sm from statsmodels.formula.api import ols
statsmodelsは初出です。
Rのフォーミュラ式ライクにモデルを作成することが可能です。
olsはOrdinary Least Squares regression:最小二乗回帰ですね。
一度モデルにしてから、anovaにかけます。
データの読み込み
本書に書いてあるとおり、
npkデータを使用します。
これはエンドウマメの栽培データで、説明変数として畑の場所と与えた肥料における窒素・リン・カリウムの組み合わせが、目的変数として1/70エーカー辺りの収穫量（ポンド)が与えられています。
とのこと
df = pd.read_csv("npl.csv")
df.head()
block N P K yield 0 1 0 1 1 49.5 1 1 1 1 0 62.8 2 1 0 0 0 46.8 3 1 1 0 1 57.0 4 2 1 0 0 59.8
anova
# カラム名がyieldだとエラーになるので変えておく df.columns = ["block","N","P","K","yield_col"] # 線形回帰のモデルをつくる npk_lm = ols("yield_col~block*N*P*K",data=df).fit() # anova result = sm.stats.anova_lm(npl_lm) result
df sum_sq mean_sq F PR(>F) block 5 343.295000 68.659000 4.446666 0.015939 N 1 189.281667 189.281667 12.258734 0.004372 P 1 8.401667 8.401667 0.544130 0.474904 K 1 95.201667 95.201667 6.165689 0.028795 N:P 1 21.281667 21.281667 1.378297 0.263165 N:K 1 33.135000 33.135000 2.145972 0.168648 P:K 1 0.481667 0.481667 0.031195 0.862752 N:P:K 1 1.651002 1.651002 0.106926 0.749304 Residual 12 185.286667 15.440556 NaN NaN
ほぼ本書と同様の結果です。
なお、stats.anova_lmは引数にtyp=タイプを与えることで、
anovaのtypeを指定可能です。
前回はこちら esu-ko.hatenablog.com
