Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード Twitter

2017-10-03

R・Pythonのロジスティック回帰結果が一致することを確認

| 17:02 | R・Pythonのロジスティック回帰結果が一致することを確認を含むブックマーク

ちゃんと一致するのかなって不安になったので適当なデータでチェック。

元になるデータはこの記事のものをそのまま流用している&何なのかは知らん。

Python

その1

Pythonのコード

import pandas as pd
import sklearn.linear_model
df = pd.read_csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
# rank列をone-hot encodingで0・1に変換&はじめの列は余計なので削除(drop_first=True)&元のrank変数もいらないのでカット
df = pd.concat([df, pd.get_dummies(df['rank'], prefix='rank', drop_first=True)], axis=1)
df.drop('rank', axis=1, inplace=True)

# 説明変数(x)と非説明変数(y)に分ける
y = df['admit'].values
x = df.loc[:, df.columns != 'admit'].values
# 学習(正則化なしにしたいならでかいCにするのか・・・ https://github.com/scikit-learn/scikit-learn/issues/6738)
clf = sklearn.linear_model.LogisticRegression(C=1000000)
clf.fit(x, y)

実行結果(回帰係数表示)

clf.coef_
Out[59]: array([[ 0.00225514,  0.80176426, -0.67456794, -1.33708169, -1.55029739]])
clf.intercept_
Out[60]: array([-3.9776519])

デフォルト正則化が入ってて、かつ、それを消すのにC=0じゃなくてC=でかい値にしなきゃいけないのがアレだなって思うんだ・・・同様に思ってる人もいるわけだ、Issueは開きっぱで放置されてるけど。。。

その2

Pythonのコード

import statsmodels.api as sm
df['intercept'] = 1.0
logit = sm.Logit(df['admit'], df.loc[:, df.columns != 'admit'])
result = logit.fit()

実行結果

result.summary()
Out[36]: 
<class 'statsmodels.iolib.summary.Summary'>
"""
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                  admit   No. Observations:                  400
Model:                          Logit   Df Residuals:                      394
Method:                           MLE   Df Model:                            5
Date:                Tue, 03 Oct 2017   Pseudo R-squ.:                 0.08292
Time:                        21:28:21   Log-Likelihood:                -229.26
converged:                       True   LL-Null:                       -249.99
                                        LLR p-value:                 7.578e-08
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
gre            0.0023      0.001      2.070      0.038       0.000       0.004
gpa            0.8040      0.332      2.423      0.015       0.154       1.454
rank_2        -0.6754      0.316     -2.134      0.033      -1.296      -0.055
rank_3        -1.3402      0.345     -3.881      0.000      -2.017      -0.663
rank_4        -1.5515      0.418     -3.713      0.000      -2.370      -0.733
intercept     -3.9900      1.140     -3.500      0.000      -6.224      -1.756
==============================================================================
"""

切片は自分でたさないとだめぽいのがあれだわ。

R編

Rのコード

df <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv", stringsAsFactors=FALSE)
# rankは因子扱い
df$rank <- as.factor(df$rank)
# ロジスティック回帰実行
glm(admit~., data=df, family=binomial)

実行結果

Call:  glm(formula = admit ~ ., family = binomial, data = df)

Coefficients:
(Intercept)          gre          gpa        rank2        rank3        rank4  
  -3.989979     0.002264     0.804038    -0.675443    -1.340204    -1.551464  

Degrees of Freedom: 399 Total (i.e. Null);  394 Residual
Null Deviance:	    500 
Residual Deviance: 458.5 	AIC: 470.5

まとめ

上の結果を読むと、係数は大体等しくなっているので一安心。どう見てもRの方がコードがシンプルなんだが、あえて言うなら

df.hist()

と書くと全部のヒストグラムを出してくれるのがPandasのいいところだなって思いました(まる)。

Pythonのもっと良い書き方を学びたい。

結局、sklearnとstatmodelsはどう使い分けるんだ・・・

f:id:teramonagi:20171003161756p:image