statsmodelsで分散分析する時"typ"を"type"とミスタイプすると起きること

この記事の内容

TechAcademyのデータサイエンスコース
を受講する中で躓いた点を備忘録として記載します。
このコースでは、Statsmodelsというライブラリを使って統計分析を行います。

コースは以下の書籍を読みながら進めていくスタイルになります。
記事の内容は第5部「正規線形モデル」の第3章「複数の説明変数を持つモデル」
に関するものです。

Pythonで学ぶあたらしい統計学の教科書

実行環境

macOS Mojave version10.14.5
AWS cloud9
Python 3.6.5
pandas 0.23.0
statsmodels 0.10.1
matplotlib 3.1.1
seaborn 0.9.0

結論

sm.stats.anova_lmを使って分散分析する際、
typ指定で"typ"を"type"とミスタイプすると、typ=1と同じ結果になる。
type=2と指定しても、エラーが出ることなくtyp=1と同じ結果になります。

F比やp値など分散分析の肝となる値が変わってくるので注意が必要です。

実際に計算してみる

実際に分散分析して、typ=1 , typ=2 , type=2 を指定した時の結果を見てみます。
データはkaggleのHouse Prices: Advanced Regression Techniquesを使います。住宅価格を予測する問題です。

まずサイトから落としてきたデータをDataFrameで読み込みます。
説明を簡便にするために以下の列だけを取り出し、図示します。

・Id　（物件にふられた番号）
・OverallQual　（全体的な質）
・YearBuilt　（建造年）
・1stFlrSF　（１F広さ）
・SalePrice　（販売価格）

import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
from matplotlib import pyplot as plt
import seaborn as sns
sns.set()
%matplotlib inline

data=pd.read_csv("train.csv") 
#↑サイトから入手した訓練データ

data=data[["Id","OverallQual","YearBuilt","1stFlrSF","SalePrice"]]
#↑Id , OverallQual ,YearBuilt ,1stFlrSF , SalePrice を取り出す

data.rename(columns={"1stFlrSF":"FlrSF1st"},inplace=True)
#↑statsmodel使用時にエラーが出るため名前を変える

data=data.sort_values(by="SalePrice")
#↑SalePrice列を小さい順に並び変える。実際の値と当てはめ値（統計モデルから推定した値）が比較し易くなる。

sns.pairplot(data)
#↑全ての組み合わせを図示

一番下の行を見ると、
SalePriceに対してOverallQual、YearBuilt、1stFlrSF　はいずれも相関があります。
ですので、いずれも統計モデルの変数に加えた方が良さそうです。
最も簡単なモデリングとして、重回帰モデルを作ります。
つまり、以下の式で実際のSalePrice値とSalePrice（当てはめ値）の差が最小になるように係数βを求めます。

SalePrice（当てはめ値）=β_0＋β_1×OverallQual＋β_2×YearBuilt＋β_3×FlrSF1st

重回帰モデルはstatsmodelsで作ることができます。

SalePrice_lm=smf.ols(formula="SalePrice~OverallQual+YearBuilt+FlrSF1st" , data=data).fit()

formula="SalePrice~OverallQual+YearBuilt+FlrSF1st"

の箇所はSalePriceを応答変数、OverallQual、YearBuilt、FlrSF1stを説明変数にする意味です。
このモデルから、当てはめ値を求めます。

prediction=SalePrice_lm.predict() #重回帰モデルから求めた当てはめ値

実際のSalePriceとSalePrice（当てはめ値）を図示します。

plt.scatter(range(len(data["SalePrice"])),data["SalePrice"] , label="real" , color="blue")
plt.scatter(range(len(data["SalePrice"])),SalePrice_lm.predict() , label="prediction" , color="lightblue" , alpha=0.3)

plt.ylim(0,700000)
plt.xlabel("Id")
plt.ylabel("SalePrice($)")
plt.legend(fontsize=15)

OverallQual、YearBuilt、FlrSF1stだけを説明変数とした重回帰モデルでもそこそこフィッティングできました。高級物件がうまくフィッティングできていないのは、そういう物件はOverallQual、YearBuilt、FlrSF1st以外で価値が決まっているからかもしれません。

本題に入ります。重回帰モデルを使って分散分析を行います。

typeIの分散分析表

sm.stats.anova_lm(SalePrice_lm , typ=1).round(3)

typeIIの分散分析表

sm.stats.anova_lm(SalePrice_lm , typ=2).round(3)

typeIIの分散分析表で"typ"を"type"に変更します。

sm.stats.anova_lm(SalePrice_lm , type=2).round(3)

エラーが出ることなくtyp=1とした場合と同じ結果が表示されました。
ちなみに、typを指定しないとtyp=1として計算されます。

sm.stats.anova_lm(SalePrice_lm).round(3)

ご覧のようにtypeIとtypeIIでは分析結果が変わってきます。
特に、F比と呼ばれる値やp値（図中PR(>F)1)が変わってくると判断を誤るので注意が必要です。
F比はざっくり言うと、「その変数のフィッティングへの寄与度」になります。
この値が大きいほど、その変数がないと上手くフィッティングできないと言えます。
変数を取捨選択する際に考慮する値になります。

おわり

今後もTechAcademyのデータサイエンスコースを進める中で、
躓いた点を記録していこうと思います。