PythonのStatsmodelsを使って回帰分析を行う
Pythonを使って回帰分析を行う。使用するライブラリはStatsmodelsである。
%matplotlib inline
まず対象となるデータを読み込む。これはR処理系に付属しているattitudeというデータを write.csv(attitude, "attitude.csv", quote=FALSE, row.names=FALSE)
でCSVにしたものである。
この記事を書くにあたって参考にしている「定性的データ分析」(金明哲)によると、 このデータは35人の事務員に対するアンケート結果で、それぞれの列は次のような意味である。 (ところどころ意味が分からないのだが)
rating: 総合評価, complaints: 従業員の苦情の取り扱い, privileges: 特権は許さない, learning: 学習の機会, raises: 能力に基づいた昇給 , critical: 加重, advance: 昇進
まずは散布図行列を描く。
from pandas.tools.plotting import scatter_matrix
scatter_matrix(attitude, figsize=(7, 7))
plt.show()
数字は満足度の割合(パーセンテージ)である。総合評価ratingとcomplaintsに相関が高そうだということが見て取れる。 対角線上にはその列(行といっても同じだが)の変数のヒストグラムが表示される。
次にStatsmodelsを使い、ratingを目的変数、残りのすべてを説明変数として重回帰分析を行う。 add_constantという関数は1という値だけの列を追加していて、これが切片(説明変数にかかわらずオフセットされる量)となる。 Statsmodelsの流儀で、モデルに切片を使う場合はこのようにしなければならない。
import statsmodels.api as sm
# Fit and summarize OLS model
y = attitude['rating']
X = attitude[['complaints', 'privileges', 'learning', 'raises', 'critical', 'advance']]
mod = sm.OLS(y, sm.add_constant(X))
res = mod.fit()
print res.summary()
この結果が妥当だとすると、回帰係数coefにより、complaintsやlearningが大きいほど総合評価ratingは上がる。逆にadvanceが高いほど総合評価は下がるということになるようだ。
回帰モデルの評価指標としてサマリーの右上に表示されているような指標値がある。
Log-LikelihoodとR-squaredはモデル作成に使ったデータへの当てはまりの良さを示す。 しかし一般的に説明変数が多いほどデータへの当てはまりは向上するが、それが未知のデータの予測に良いとは限らず、これらの値が高くても無駄な変数を含むかもしれない。 AICやAdj. R-squaredは当てはまりの良さを説明変数の数で打ち消すような効果を持つ。
回帰変数に対して表示されているP>|t|は回帰変数に対してt検定を行った場合のp値であり、小さいほどその変数の影響が偶然でない可能性が高い (検定を行う場合は0.05未満や0.01未満が基準となる)。 そこで、P>|t|の小さい順に変数を追加していって、評価指標がどのように変化するかを見てみる。 ここでは変数を追加していくコードを自分で書いているが、おそらく自動で行うような関数もあるのかもしれない。
def olc(i):
y = attitude['rating']
X_labels = ['complaints', 'learning', 'advance', 'privileges', 'raises', 'critical'][0:i+1]
X = attitude[X_labels]
mod = sm.OLS(y, sm.add_constant(X))
res = mod.fit()
return [res.llf, res.aic, res.rsquared, res.rsquared_adj, X_labels]
DataFrame(
[olc(i) for i in range(0, 6)],
columns=['Log-Likelihood', 'AIC', 'R-squared', 'Adj. R-squared', 'vars'])
Log-LikelihoodとAdj. R-squared(いずれも高いほうが良い)は変数を増やすのにつれて単調増加している。
一方AIC(低いほど良い)とAdj. R-squared(高いほうが良い)は向上が頭打ちになるポイントがある。 AICの場合[complaints, learning]を説明変数に使ったときが最もよく、Adj. R-squaredの場合[complaints, learning, advance]を使ったときが最もよい。
なおAICの値は「定性的データ分析」のRの実行結果に出てくる値とは一致しておらず、何か計算方法に違いがあるのかもしれない。
AICをもとにすると[complaints, learning]だけを説明変数にするのがよいようだ。改めてサマリーを生成したものを次に示す。
# Fit and summarize OLS model
y = attitude['rating']
X = attitude[['complaints', 'learning']]
mod = sm.OLS(y, sm.add_constant(X))
res = mod.fit()
print res.summary()
import statsmodels.formula.api as smf
mod = smf.ols(formula="rating ~ complaints + learning", data=attitude)
res = mod.fit()
print res.summary()
回帰分析では説明変数同士の強い相関によって回帰モデルに矛盾が生じることがある。 これをチェックする指標としてVIF (Variance Inflation Factor) という値がある。
StatsmodelsでVIFを計算する関数variance_inflation_factorは少し使いにくい形をしている。 説明変数とVIFを計算したい列の番号を引数として与える。 ここでは切片を含めてすべての列に対するVIFを計算した。
num_cols = mod.exog.shape[1] # 説明変数の列数
vifs = [variance_inflation_factor(mod.exog, i)
for i in range(0, num_cols)]
DataFrame(vifs, index=mod.exog_names, columns=["VIF"])
VIFの値が10より大きいときは多重共線性があるとされる。 Statsmodels の variance_inflation_factor のリファレンス では5より大きい場合となっている。 この場合 complaints も learning も1.56程度であり、多重共線性はないと判断される。切片の分はよくわからないが、切片についてVIFを計算している例をみたことがないので無視してよいのだろう。