pythonで一般化線形モデルと公式リファレンスの罠 - using statsmodels -
寒くなってきました。最近、pythonでデータの解析をすることにいそしんでおります。
Rでできることをpythonでやりたいなと思っていろいろ調べてみると、まぁなかなかできるようになっていなかったりするわけで、その辺を整備し始めたので、ここに書いていこうと思います。
一般化線形モデル（Generalized Linear Model）
一般化線形モデルは、回帰の問題ではよく利用される方法で、Rではglmという関数で使うことができます。有名なものはロジスティック回帰です。図的には、以下のような感じになります。
これをpythonで行うときには、多くの場合scikit-learnを使うような気がするのですが、scikit-learnのアウトプットは、統計的に結果を見たい場合に十分とはいえません。できればRのロジスティック回帰の結果みたいなの（下図）、、、を出力したいんです。
statsmodelsを使ってみよう。
そこで、そんな要望に答えるために、statsmodelsというモジュールが提供されています。どうもこれを使用すれば、Rのglm的なコトができるらしいと聞きつけて、やってみました。以下のコードをのせているので興味のある人は是非やってみて欲しいと思います。最後にRの対応するコードも付けておきました。使用しているデータはstatsmodelsの中のものですが、以下のリンクからもダウンロードできるようにしました。
https://dl.dropboxusercontent.com/u/2133906/sample.csv
※最後に注意事項も書いていますので、そこもできれば、読んでください。
## 一般化線形モデル
from __future__ import print_function
import numpy as np
import pandas as pd
from pandas import DataFrame, Series
from matplotlib import pyplot as plt
#glmを可能にするために必要なものなど
#statmodelsのGLMを使えるようにするモジュール
import statsmodels.api as sm
#GLMの中で用いる統計処理関数のインポート, chi**2の値などの算出に使用している
from scipy import stats
#Rのglmを使用可能にするための関数
import statsmodels.formula.api as smf
# ## OLS, GLM: Gaussian response data
#sampleデータのロード
#形式はpandas.DataFrame, naは削除する必要あり
df = sm.datasets.get_rdataset("Guerry", "HistData").data
df = df[['Lottery', 'Literacy', 'Wealth', 'Region']].dropna()
df.head()
#線形回帰モデルの構築, 最小自乗法を使う場合
#modelの構築
ols_model = 'Lottery ~ Literacy + Wealth + Region'
#モデルの構築
mod = smf.ols(formula=ols_model, data=df)
#fitに回帰した結果が入っているので、これをresに代入する
res = mod.fit()
#結果の参照
res.summary()
"""
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: Lottery No. Observations: 85
Model: GLM Df Residuals: 78
Model Family: Gaussian Df Model: 6
Link Function: identity Scale: 436.433771252
Method: IRLS Log-Likelihood: -375.30
Date: Fri, 07 Nov 2014 Deviance: 34042.
Time: 18:59:36 Pearson chi2: 3.40e+04
No. Iterations: 4
===============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
-------------------------------------------------------------------------------
Intercept 38.6517 9.456 4.087 0.000 20.118 57.186
Region[T.E] -15.4278 9.727 -1.586 0.113 -34.493 3.637
Region[T.N] -10.0170 9.260 -1.082 0.279 -28.167 8.133
Region[T.S] -4.5483 7.279 -0.625 0.532 -18.815 9.718
Region[T.W] -10.0913 7.196 -1.402 0.161 -24.195 4.013
Literacy -0.1858 0.210 -0.886 0.376 -0.597 0.225
Wealth 0.4515 0.103 4.390 0.000 0.250 0.653
===============================================================================
"""
#線形回帰モデルの構築, glmによるニュートン法を用いる場合
#modelの構築
glm_gauss_model = 'Lottery ~ Literacy + Wealth + Region'
#モデルの構築
mod = smf.glm(formula=glm_gauss_model, data=df, family=sm.families.Gaussian())
#fitに回帰した結果が入っているので、これをresに代入する
res = mod.fit()
#結果の参照
res.summary()
# ## GLM: Binomial response data
#
# ### Load data
#
# In this example, we use the Star98 dataset which was taken with permission
# from Jeff Gill (2000) Generalized linear models: A unified approach. Codebook
# information can be obtained by typing:
#一応データセットに関するディスクリプションを参照しておく。
print(sm.datasets.star98.NOTE)
"""
簡単に言えば、目的変数に「数学の成績が国の平均以上」だった人と、
「以下だった人」の数がたぶん学校毎に入っていて、それを学校に関係する
その他の変数で説明しましょうってことだと思います。
"""
#データセットの読み込み
star98 = sm.datasets.star98.load()
#pandas形式にしたい場合は
#star98 = sm.datasets.star98.load_pandas().data
"""
大事：
本来なら上記のロードだけで良いのだが、、、
pythonのstatsmodelsの通常のGLM = sm.GLMは
どういうわけか、intercept項を入れないまま回帰するので
sm.GLMを用いる場合には定数項を入れる必要があります。
"""
#data.exog = sm.add_constant(data.exog, prepend=False)
#datasetの形式を確認
type(star98.data)
# >>> numpy.core.records.recarray
#データセットをpandasのDataFrameにしておく。
star98 = DataFrame(star98.data)
#データセットの最初の5件を確認
star98.head()
#一般化線形モデルの構築
#2項回帰モデル、link:logit
"""
Binomialでモデリングする際には、
通常のRでは目的変数にcbind(成功数, 失敗数)を取るようにする。
これをpythonのstatsmodelsでは、~の前に以下のように記述する
'成功数 + 失敗数 ~ 説明変数'
"""
#モデルを定義する。
formula_binom = 'NABOVE + NBELOW ~ LOWINC+ PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF'
#成功割合の変数を作成する
star98['SUCCESS'] = star98['NABOVE'] / (star98['NABOVE'] + star98['NBELOW'])
star98['TOTAL'] = star98['NABOVE'] + star98['NBELOW']
# binomialの回帰モデルを構築する
glm_binom = smf.glm(formula=formula_binom, data=star98, family=sm.families.Binomial())
# res_binomにglmの結果を保存する
res_binom = glm_binom.fit()
# summaryを表示
res_binom.summary()
Generalized Linear Model Regression Results
================================================================================
Dep. Variable: ['NABOVE', 'NBELOW'] No. Observations: 303
Model: GLM Df Residuals: 282
Model Family: Binomial Df Model: 20
Link Function: logit Scale: 1.0
Method: IRLS Log-Likelihood: -2998.6
Date: Fri, 07 Nov 2014 Deviance: 4078.8
Time: 19:58:26 Pearson chi2: 4.05e+03
No. Iterations: 7
============================================================================================
coef std err z P>|z| [95.0% Conf. Int.]
--------------------------------------------------------------------------------------------
Intercept 2.9589 1.547 1.913 0.056 -0.073 5.990
LOWINC -0.0168 0.000 -38.749 0.000 -0.018 -0.016
PERASIAN 0.0099 0.001 16.505 0.000 0.009 0.011
PERBLACK -0.0187 0.001 -25.182 0.000 -0.020 -0.017
PERHISP -0.0142 0.000 -32.818 0.000 -0.015 -0.013
PCTCHRT 0.0049 0.001 3.921 0.000 0.002 0.007
PCTYRRND -0.0036 0.000 -15.878 0.000 -0.004 -0.003
PERMINTE 0.2545 0.030 8.498 0.000 0.196 0.313
AVYRSEXP 0.2407 0.057 4.212 0.000 0.129 0.353
PERMINTE:AVYRSEXP -0.0141 0.002 -7.391 0.000 -0.018 -0.010
AVSALK 0.0804 0.014 5.775 0.000 0.053 0.108
PERMINTE:AVSALK -0.0040 0.000 -8.450 0.000 -0.005 -0.003
AVYRSEXP:AVSALK -0.0039 0.001 -4.059 0.000 -0.006 -0.002
PERMINTE:AVYRSEXP:AVSALK 0.0002 2.99e-05 7.428 0.000 0.000 0.000
PERSPENK -1.9522 0.317 -6.162 0.000 -2.573 -1.331
PTRATIO -0.3341 0.061 -5.453 0.000 -0.454 -0.214
PERSPENK:PTRATIO 0.0917 0.015 6.321 0.000 0.063 0.120
PCTAF -0.1690 0.033 -5.169 0.000 -0.233 -0.105
PERSPENK:PCTAF 0.0490 0.007 6.574 0.000 0.034 0.064
PTRATIO:PCTAF 0.0080 0.001 5.362 0.000 0.005 0.011
PERSPENK:PTRATIO:PCTAF -0.0022 0.000 -6.445 0.000 -0.003 -0.002
============================================================================================
# ### 興味のある統計量を表示する方法
#試行回数
print("Total number of trials:" + str(star98[['NABOVE','NBELOW']].sum().sum(axis=1)))
#パラメータの表示
print('Parameters: ', res_binom.params)
#各パラメータに対するt値
print('T-values: ', res_binom.tvalues)
# ### Plots
# ロジスティック回帰の結果について図で確認する。
# プロットを通して、データの予測精度について分析する
# ・予測値と実測値のプロット
# ・予測値と残差のプロット
#観測値
y = star98['SUCCESS']
#推定値
yhat = Series(res_binom.mu)
#linear_predictor
y_lp = Series(np.log(res_binom.mu/(1-res_binom.mu)))
#残差
residuals = Series(res_binom.resid_pearson)
#dataframe:
dd = DataFrame({
'OBS':y,
'FITTED':yhat,
'LINEAR_PRED':y_lp,
'RESIDUALS':residuals
})
# 予測値と観測値の差を確認する
from ggplot import *
p_fitted_vs_obs = ggplot(aes(x='FITTED',y='OBS'),data=dd)+geom_point(colour="red")+geom_abline(intercept=0,slope=1)
# 予測値と残差の関係性を確認する
p_fitted_vs_residuals = ggplot(aes(x='FITTED',y='RESIDUALS'),data=dd)+geom_point(colour="red")+stat_smooth()
Rのコード
以下のようになります。よって、結果が一致しました！
> star98 = read.csv("sample.csv")
>
> glm_binom = glm(cbind(NABOVE,NBELOW) ~ LOWINC + PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF, data=star98, family="binomial")
> summary(glm_binom)
Call:
glm(formula = cbind(NABOVE, NBELOW) ~ LOWINC + PERASIAN + PERBLACK +
PERHISP + PCTCHRT + PCTYRRND + PERMINTE * AVYRSEXP * AVSALK +
PERSPENK * PTRATIO * PCTAF, family = "binomial", data = star98)
Deviance Residuals:
Min 1Q Median 3Q Max
-9.5217 -2.6505 -0.3357 2.1658 15.4200
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.959e+00 1.547e+00 1.913 0.0557 .
LOWINC -1.682e-02 4.340e-04 -38.749 < 2e-16 ***
PERASIAN 9.925e-03 6.014e-04 16.505 < 2e-16 ***
PERBLACK -1.872e-02 7.436e-04 -25.182 < 2e-16 ***
PERHISP -1.424e-02 4.339e-04 -32.818 < 2e-16 ***
PCTCHRT 4.917e-03 1.254e-03 3.921 8.81e-05 ***
PCTYRRND -3.580e-03 2.255e-04 -15.878 < 2e-16 ***
PERMINTE 2.545e-01 2.995e-02 8.498 < 2e-16 ***
AVYRSEXP 2.407e-01 5.714e-02 4.212 2.53e-05 ***
AVSALK 8.041e-02 1.392e-02 5.775 7.69e-09 ***
PERSPENK -1.952e+00 3.168e-01 -6.162 7.19e-10 ***
PTRATIO -3.341e-01 6.126e-02 -5.453 4.95e-08 ***
PCTAF -1.690e-01 3.270e-02 -5.169 2.36e-07 ***
PERMINTE:AVYRSEXP -1.408e-02 1.905e-03 -7.391 1.46e-13 ***
PERMINTE:AVSALK -4.005e-03 4.740e-04 -8.450 < 2e-16 ***
AVYRSEXP:AVSALK -3.906e-03 9.624e-04 -4.059 4.92e-05 ***
PERSPENK:PTRATIO 9.171e-02 1.451e-02 6.321 2.60e-10 ***
PERSPENK:PCTAF 4.899e-02 7.452e-03 6.574 4.89e-11 ***
PTRATIO:PCTAF 8.041e-03 1.499e-03 5.362 8.22e-08 ***
PERMINTE:AVYRSEXP:AVSALK 2.220e-04 2.989e-05 7.428 1.10e-13 ***
PERSPENK:PTRATIO:PCTAF -2.249e-03 3.490e-04 -6.445 1.15e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 34345.4 on 302 degrees of freedom
Residual deviance: 4078.8 on 282 degrees of freedom
AIC: 6039.2
Number of Fisher Scoring iterations: 4
公式リファレンスの罠。
githubにあるロジスティック回帰のコードが入っていますが、これはおそらく間違いです。
statsmodels/glm_formula.py at master · statsmodels/statsmodels · GitHub
どういうわけか、以下のコマンドでロジスティック回帰が実行できますということになっていますが...
## Generalized Linear Models (Formula)
# This notebook illustrates how you can use R-style formulas to fit Generalized Linear Models.
#
# To begin, we load the ``Star98`` dataset and we construct a formula and pre-process the data:
from __future__ import print_function
import statsmodels.api as sm
import statsmodels.formula.api as smf
star98 = sm.datasets.star98.load_pandas().data
formula = 'SUCCESS ~ LOWINC + PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF'
dta = star98[['NABOVE', 'NBELOW', 'LOWINC', 'PERASIAN', 'PERBLACK', 'PERHISP',
'PCTCHRT', 'PCTYRRND', 'PERMINTE', 'AVYRSEXP', 'AVSALK',
'PERSPENK', 'PTRATIO', 'PCTAF']]
endog = dta['NABOVE'] / (dta['NABOVE'] + dta.pop('NBELOW'))
del dta['NABOVE']
dta['SUCCESS'] = endog
# Then, we fit the GLM model:
mod1 = smf.glm(formula=formula, data=dta, family=sm.families.Binomial()).fit()
mod1.summary()問題があるのは、formulaの部分です。この形式では「目的変数」にSUCCESS（成功割合）のみを取ればよいということになっていますが、ロジスティック回帰の結果というのは「どれだけの回数試行を行ったか」に依存するため、成功の割合のみを回帰すれば良いというのは間違いです。実際、この結果を確認すると、以下のようになります。
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: SUCCESS No. Observations: 303
Model: GLM Df Residuals: 282
Model Family: Binomial Df Model: 20
Link Function: logit Scale: 1.0
Method: IRLS Log-Likelihood: -189.70
Date: Fri, 07 Nov 2014 Deviance: 380.66
Time: 20:11:04 Pearson chi2: 8.48
No. Iterations: 7
============================================================================================
coef std err z P>|z| [95.0% Conf. Int.]
--------------------------------------------------------------------------------------------
Intercept 0.4037 25.036 0.016 0.987 -48.665 49.472
LOWINC -0.0204 0.010 -1.982 0.048 -0.041 -0.000
PERASIAN 0.0159 0.017 0.910 0.363 -0.018 0.050
PERBLACK -0.0198 0.020 -1.004 0.316 -0.058 0.019
PERHISP -0.0096 0.010 -0.951 0.341 -0.029 0.010
PCTCHRT -0.0022 0.022 -0.103 0.918 -0.045 0.040
PCTYRRND -0.0022 0.006 -0.348 0.728 -0.014 0.010
PERMINTE 0.1068 0.787 0.136 0.892 -1.436 1.650
AVYRSEXP -0.0411 1.176 -0.035 0.972 -2.346 2.264
PERMINTE:AVYRSEXP -0.0031 0.054 -0.057 0.954 -0.108 0.102
AVSALK 0.0131 0.295 0.044 0.965 -0.566 0.592
PERMINTE:AVSALK -0.0019 0.013 -0.145 0.885 -0.028 0.024
AVYRSEXP:AVSALK 0.0008 0.020 0.038 0.970 -0.039 0.041
PERMINTE:AVYRSEXP:AVSALK 5.978e-05 0.001 0.068 0.946 -0.002 0.002
PERSPENK -0.3097 4.233 -0.073 0.942 -8.606 7.987
PTRATIO 0.0096 0.919 0.010 0.992 -1.792 1.811
PERSPENK:PTRATIO 0.0066 0.206 0.032 0.974 -0.397 0.410
PCTAF -0.0143 0.474 -0.030 0.976 -0.944 0.916
PERSPENK:PCTAF 0.0105 0.098 0.107 0.915 -0.182 0.203
PTRATIO:PCTAF -0.0001 0.022 -0.005 0.996 -0.044 0.044
PERSPENK:PTRATIO:PCTAF -0.0002 0.005 -0.051 0.959 -0.010 0.009
============================================================================================
回帰係数が異なっている分には実際、問題はなくて、問題があるのは「Yの推定値」に差があるかどうかです。それを確認するために、結果を確認します。まず、正しい結果はこのようになっています。
In [280]: yhat.head() Out[280]: 0 0.583312 1 0.751447 2 0.500583 3 0.685345 4 0.322510 dtype: float64
一方で間違っている結果は、このようになります。
In [279]: Series(mod1.mu).head() Out[279]: 0 0.578461 1 0.769802 2 0.449346 3 0.683891 4 0.247756 dtype: float64
先ほども説明した通り、公式リファレンスの罠ですね。下の結果は正しいロジスティック回帰の結果ではありません。この原因を作っているのは、先ほども言ったように成功割合のみでロジスティック回帰を実行したせいです。なぜなら、本来の成功回数Y(i)は、試行数N(i)、成功確率Pi(i)の2項分布に従います。しかしながら、Y(i)/N(i)の分布は当たり前ですが2項分布ではありません。確率変数変換をすれば自明です。なぜなら、P(i)の値域は0-1ですからね。結果的に、ロジスティック回帰は結果が違うことになってしまうというわけです。
注意）実際、成功割合のみを与えるようなロジスティック回帰が実行される場合には「全てのケースでの試行回数」が同じであるという仮定のもとで実行されていますので、結果は異なって当然です。
なので、pythonで行う際にはその点はご注意ください。
最後に
今回、pythonでglmを行う方法を考えましたが、まぁ実際Rでしようよっていうのが正直なところです。残念なことに笑
