PythonのStatsmodelsを使って回帰分析を行う

 

Pythonを使って回帰分析を行う。使用するライブラリはStatsmodelsである。

In [78]:
%matplotlib inline
 

まず対象となるデータを読み込む。これはR処理系に付属しているattitudeというデータを write.csv(attitude, "attitude.csv", quote=FALSE, row.names=FALSE)CSVにしたものである。

In [79]:
attitude = pd.read_csv('attitude.csv')
attitude.columns
Out[79]:
Index([u'rating', u'complaints', u'privileges', u'learning', u'raises',
       u'critical', u'advance'],
      dtype='object')
 

この記事を書くにあたって参考にしている「定性的データ分析」(金明哲)によると、 このデータは35人の事務員に対するアンケート結果で、それぞれの列は次のような意味である。 (ところどころ意味が分からないのだが)

rating: 総合評価, complaints: 従業員の苦情の取り扱い, privileges: 特権は許さない, learning: 学習の機会, raises: 能力に基づいた昇給 , critical: 加重, advance: 昇進

まずは散布図行列を描く。

In [80]:
from pandas.tools.plotting import scatter_matrix
scatter_matrix(attitude, figsize=(7, 7))
plt.show()
 
 

数字は満足度の割合(パーセンテージ)である。総合評価ratingとcomplaintsに相関が高そうだということが見て取れる。 対角線上にはその列(行といっても同じだが)の変数のヒストグラムが表示される。

次にStatsmodelsを使い、ratingを目的変数、残りのすべてを説明変数として重回帰分析を行う。 add_constantという関数は1という値だけの列を追加していて、これが切片(説明変数にかかわらずオフセットされる量)となる。 Statsmodelsの流儀で、モデルに切片を使う場合はこのようにしなければならない。

In [81]:
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()
 
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 rating   R-squared:                       0.733
Model:                            OLS   Adj. R-squared:                  0.663
Method:                 Least Squares   F-statistic:                     10.50
Date:                Tue, 08 Nov 2016   Prob (F-statistic):           1.24e-05
Time:                        23:12:33   Log-Likelihood:                -97.250
No. Observations:                  30   AIC:                             208.5
Df Residuals:                      23   BIC:                             218.3
Df Model:                           6                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         10.7871     11.589      0.931      0.362       -13.187    34.761
complaints     0.6132      0.161      3.809      0.001         0.280     0.946
privileges    -0.0731      0.136     -0.538      0.596        -0.354     0.208
learning       0.3203      0.169      1.901      0.070        -0.028     0.669
raises         0.0817      0.221      0.369      0.715        -0.376     0.540
critical       0.0384      0.147      0.261      0.796        -0.266     0.342
advance       -0.2171      0.178     -1.218      0.236        -0.586     0.152
==============================================================================
Omnibus:                        2.386   Durbin-Watson:                   1.795
Prob(Omnibus):                  0.303   Jarque-Bera (JB):                1.255
Skew:                          -0.081   Prob(JB):                        0.534
Kurtosis:                       2.011   Cond. No.                     1.34e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.34e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
 

この結果が妥当だとすると、回帰係数coefにより、complaintsやlearningが大きいほど総合評価ratingは上がる。逆にadvanceが高いほど総合評価は下がるということになるようだ。

回帰モデルの評価指標としてサマリーの右上に表示されているような指標値がある。

Log-LikelihoodとR-squaredはモデル作成に使ったデータへの当てはまりの良さを示す。 しかし一般的に説明変数が多いほどデータへの当てはまりは向上するが、それが未知のデータの予測に良いとは限らず、これらの値が高くても無駄な変数を含むかもしれない。 AICやAdj. R-squaredは当てはまりの良さを説明変数の数で打ち消すような効果を持つ。

回帰変数に対して表示されているP>|t|は回帰変数に対してt検定を行った場合のp値であり、小さいほどその変数の影響が偶然でない可能性が高い (検定を行う場合は0.05未満や0.01未満が基準となる)。 そこで、P>|t|の小さい順に変数を追加していって、評価指標がどのように変化するかを見てみる。 ここでは変数を追加していくコードを自分で書いているが、おそらく自動で行うような関数もあるのかもしれない。

In [82]:
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'])
Out[82]:
 Log-LikelihoodAICR-squaredAdj. R-squaredvars
0 -99.881922 203.763844 0.681314 0.669933 [complaints]
1 -98.569362 203.138724 0.708015 0.686387 [complaints, learning]
2 -97.637912 203.275825 0.725595 0.693933 [complaints, learning, advance]
3 -97.431718 204.863435 0.729341 0.686036 [complaints, learning, advance, privileges]
4 -97.294306 206.588612 0.731809 0.675936 [complaints, learning, advance, privileges, ra...
5 -97.249909 208.499817 0.732602 0.662846 [complaints, learning, advance, privileges, ra...
 

Log-LikelihoodとAdj. R-squared(いずれも高いほうが良い)は変数を増やすのにつれて単調増加している。

一方AIC(低いほど良い)とAdj. R-squared(高いほうが良い)は向上が頭打ちになるポイントがある。 AICの場合[complaints, learning]を説明変数に使ったときが最もよく、Adj. R-squaredの場合[complaints, learning, advance]を使ったときが最もよい。

なおAICの値は「定性的データ分析」のRの実行結果に出てくる値とは一致しておらず、何か計算方法に違いがあるのかもしれない。

AICをもとにすると[complaints, learning]だけを説明変数にするのがよいようだ。改めてサマリーを生成したものを次に示す。

In [83]:
# 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()
 
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 rating   R-squared:                       0.708
Model:                            OLS   Adj. R-squared:                  0.686
Method:                 Least Squares   F-statistic:                     32.74
Date:                Tue, 08 Nov 2016   Prob (F-statistic):           6.06e-08
Time:                        23:12:33   Log-Likelihood:                -98.569
No. Observations:                  30   AIC:                             203.1
Df Residuals:                      27   BIC:                             207.3
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const          9.8709      7.061      1.398      0.174        -4.618    24.359
complaints     0.6435      0.118      5.432      0.000         0.400     0.887
learning       0.2112      0.134      1.571      0.128        -0.065     0.487
==============================================================================
Omnibus:                        6.448   Durbin-Watson:                   1.958
Prob(Omnibus):                  0.040   Jarque-Bera (JB):                1.959
Skew:                          -0.041   Prob(JB):                        0.375
Kurtosis:                       1.751   Cond. No.                         503.
==============================================================================

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

Statsmodelsにはstatsmodels.formula.apiというモジュールもある。 これを使うと目的変数と説明変数の間の関係を式で表現することができる。 この式はRのlm関数で使われるものと同じなので、Rで書かれた本などをPythonに翻訳するときに便利だ。

In [84]:
import statsmodels.formula.api as smf

mod = smf.ols(formula="rating ~ complaints + learning", data=attitude)
res = mod.fit()
print res.summary()
 
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 rating   R-squared:                       0.708
Model:                            OLS   Adj. R-squared:                  0.686
Method:                 Least Squares   F-statistic:                     32.74
Date:                Tue, 08 Nov 2016   Prob (F-statistic):           6.06e-08
Time:                        23:12:33   Log-Likelihood:                -98.569
No. Observations:                  30   AIC:                             203.1
Df Residuals:                      27   BIC:                             207.3
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept      9.8709      7.061      1.398      0.174        -4.618    24.359
complaints     0.6435      0.118      5.432      0.000         0.400     0.887
learning       0.2112      0.134      1.571      0.128        -0.065     0.487
==============================================================================
Omnibus:                        6.448   Durbin-Watson:                   1.958
Prob(Omnibus):                  0.040   Jarque-Bera (JB):                1.959
Skew:                          -0.041   Prob(JB):                        0.375
Kurtosis:                       1.751   Cond. No.                         503.
==============================================================================

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

回帰分析では説明変数同士の強い相関によって回帰モデルに矛盾が生じることがある。 これをチェックする指標としてVIF (Variance Inflation Factor) という値がある。

StatsmodelsでVIFを計算する関数variance_inflation_factorは少し使いにくい形をしている。 説明変数とVIFを計算したい列の番号を引数として与える。 ここでは切片を含めてすべての列に対するVIFを計算した。

In [85]:
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"])
Out[85]:
 VIF
Intercept 32.190129
complaints 1.553021
learning 1.553021
 

VIFの値が10より大きいときは多重共線性があるとされる。 Statsmodels の variance_inflation_factor のリファレンス では5より大きい場合となっている。 この場合 complaints も learning も1.56程度であり、多重共線性はないと判断される。切片の分はよくわからないが、切片についてVIFを計算している例をみたことがないので無視してよいのだろう。

In [ ]: