Python "statsmodels" は全体的に stable で線形回帰計算 OLS でもお世話になっているが，細かく見ていくと？？？ となるケースも発生する．ここでは，ロジステック回帰に関する２つのやり方を挙げ，注意点をメモする．
対象は **statsmodels (0.6.1) **で現時点でのLatest Stable Releaseである．
Logit Model を使うやり方
googleで検索すると，stackoverflow.com Q&Aで出てくるのが，Logit Modelを使うやり方である．Statsmodels ドキュメントでも"Regression with Discrete Dependent Variable" のところにきちんと書かれているので，ドキュメントに従って動きを確認した．
例題として解いているのは Spector and Mazzeo のデータ分析で，出典は W. Greene. "Econometric Analysis" Prentice Hall, 5th. edition. 2003 ，計量経済学の教科書である．説明変数：psi, tuce, gpa から被説明変数：post gradeを予測する問題．インターネットから得た情報から分かりやすく（独断で）置き換えると，１学期の通知表(gpa)，夏季全国模試の成績(tuce)，夏期講習への参加状況(psi)から２学期の成績が上がる/下がる(post grade)を推定するという問題である．Statsmodels ドキュメントに従って以下のようにした．
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
# Load the data from Spector and Mazzeo (1980)
spector_data = sm.datasets.spector.load()
spector_data.exog = sm.add_constant(spector_data.exog)
# Follow statsmodles ipython notebook
logit_mod = sm.Logit(spector_data.endog, spector_data.exog)
logit_res = logit_mod.fit(disp=0)
print('Parameters: ', logit_res.params)
print logit_res.summary()
上記の通り，sm.Logit() を使って問題なく計算できた．summary() は以下の通りとなった．
Logit Regression Results
==============================================================================
Dep. Variable: y No. Observations: 32
Model: Logit Df Residuals: 28
Method: MLE Df Model: 3
Date: Tue, 01 Sep 2015 Pseudo R-squ.: 0.3740
Time: 22:20:41 Log-Likelihood: -12.890
converged: True LL-Null: -20.592
LLR p-value: 0.001502
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const -13.0213 4.931 -2.641 0.008 -22.687 -3.356
x1 2.8261 1.263 2.238 0.025 0.351 5.301
x2 0.0952 0.142 0.672 0.501 -0.182 0.373
x3 2.3787 1.065 2.234 0.025 0.292 4.465
==============================================================================
GLM (Generalized Linear Models) を使うやり方
ロジスティック回帰は，リンク関数にLogit() を用いる一般化回帰モデル(GLM)の一種であるので，以下のようにも書くことができる．
# logit_mod = sm.Logit(spector_data.endog, spector_data.exog)
glm_model = sm.GLM(spector_data.endog, spector_data.exog, family=sm.families.Binomial())
# logit_res = logit_mod.fit(disp=0)
glm_reslt = glm_model.fit()
# print logit_res.summary()
print glm_reslt.summary()
sm.GLM() で family オプションをつけて使用する分布がBinomial となる指示をする．これで，リンク関数は(Binomialのディフォルトの) logit() を使う処理となる．出力される summary() は次の通り．
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 32
Model: GLM Df Residuals: 28
Model Family: Binomial Df Model: 3
Link Function: logit Scale: 1.0
Method: IRLS Log-Likelihood: -12.890
Date: Tue, 01 Sep 2015 Deviance: 25.779
Time: 22:21:20 Pearson chi2: 27.3
No. Iterations: 7
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const -13.0213 4.931 -2.641 0.008 -22.687 -3.356
x1 2.8261 1.263 2.238 0.025 0.351 5.301
x2 0.0952 0.142 0.672 0.501 -0.182 0.373
x3 2.3787 1.065 2.234 0.025 0.292 4.465
==============================================================================
Logit()を使ったケースと比較すると出力のフォーマットが異なっている．当初は，Logit() と GLM() は，wrapper が異なっているだけで(aliasingされているのであって），中の処理は一緒と考えていた．しかし，出力のMethodのところをよく見ると，
- Logit(): MLE ... most likelihood esitimation
- GLM(): IRLS ... iteratively reweighted least squares
のように，異なっている．これは別物か？とwikipediaの記事を参照しながらかなり悩んだが，
- IRLSもLikelihood（尤度）を計算する手法の一つであること．
- 算定されたパラメータや統計情報量の数値が両者で全く同じであること．
から，両者は wrapper が異なっているのみで中の処理は同じではないかと推定される．
２つの方法の比較 - 残差計算
回帰分析の状況を把握するために残差を計算することにした．OLSのやり方に従いやってみると，
>>> resid1 = logit_res.resid
AttributeError: 'LogitResults' object has no attribute 'resid'
>>> resid2 = glm_reslt.resid
AttributeError: 'GLMResults' object has no attribute 'resid'
どちらもAttributeErrorを返される．「残差」のコンセプトが線形モデルと異なるので，同じattribute名にしなかったのかも知れない．オブジェクトのクラス名から必死でドキュメントを調べると次のような正解にたどり着いた．
>>> resid1 = logit_res.resid_dev # for Logit model
>>> resid1
array([-0.23211021, -0.35027122, -0.64396264, -0.22909819, 1.06047795,
-0.26638437, -0.23178275, -0.32537884, -0.48538752, 0.85555565,
-0.22259715, -0.64918082, -0.88199929, 1.81326864, -0.94639849,
-0.24758297, -0.3320177 , -0.28054444, -1.33513084, 0.91030269,
-0.35592175, 0.44718924, -0.74400503, -1.95507406, 0.59395382,
1.20963752, 0.95233204, -0.85678568, 0.58707192, 0.33529199,
-1.22731092, 2.09663887])
>>> resid2 = glm_reslt.resid_deviance # for GLM model
>>> resid2
array([-0.23211021, -0.35027122, -0.64396264, -0.22909819, 1.06047795,
-0.26638437, -0.23178275, -0.32537884, -0.48538752, 0.85555565,
-0.22259715, -0.64918082, -0.88199929, 1.81326864, -0.94639849,
-0.24758297, -0.3320177 , -0.28054444, -1.33513084, 0.91030269,
-0.35592175, 0.44718924, -0.74400503, -1.95507406, 0.59395382,
1.20963752, 0.95233204, -0.85678568, 0.58707192, 0.33529199,
-1.22731092, 2.09663887])
両者は同じ内容なのに classによってattributeの名前が異なることについては，かなり違和感がある．（書いているProgrammerが違うのかもしれません．）
但し，出力されたresidualsの数値は全く同一で，これは回帰計算の処理が同じであったことを裏付けていると考えられる．
**Statsmodels (0.6.1)**もまだまだbeta versionなのかも知れないが，将来versionではもう少し整理して，分かりやすい体系にして欲しいと思う．
参考文献 (web site)
- Statsmodels domumentation http://statsmodels.sourceforge.net/stable/index.html
- Wikipedia "Iteratively reweighted least squares(IRLS)"
https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares