読者です 読者をやめる 読者になる 読者になる

old school magic

機械学習と統計とプログラミングについてちょっとずつ勉強していきます。

Pythonで一般化線形モデル

Python 統計 機械学習

概要

統計の勉強の一環で、最近はこの本を読んでます。かなり分かりやすいです。

統計モデリングに関する本です。一般化線形モデルを中心に話が進んでいきます。
この本はRを中心に話が進んでいきますが、せっかくなのでPythonで一般化線形モデルを試してみようと思います。

Pythonの統計ライブラリ

一般化線形モデルとは、線形回帰を(正規分布以外でも使えるように)拡張した統計モデルです。
詳しい説明は教科書に譲るとして、Pythonでのライブラリについてお話します。

Pythonで一般化線形モデル、というか統計全般のライブラリとして、簡単に調べたところメジャーなのが二種類あります。
scikit-learnStatsModels です。

scikit-learn
http://scikit-learn.org/stable/index.html

StatsModels
http://statsmodels.sourceforge.net/devel/index.html

scikit-learn はどちらかというと機械学習の、StatsModels はどちらかというと統計のライブラリです。
今回はStatsModels を使ってみようと思います。

インストール

Anaconda入れるともれなくインストールされます。
Anaconda(というかPython)のインストールはこちらを。
MacでPythonの機械学習環境構築(2014年5月版) - old school magic

Pandasというデータ分析ライブラリも(データを読み込む時に)使います。
これもAnacondaに同梱されています。

使うデータ

先程紹介した本の著者である久保先生のデータをお借りします。
http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv

どういったデータか簡単に説明すると、ある植物が種子を飛ばす数についてのデータです。
観測変数は3つあり、

  • y : 飛ばした種子の数
  • x : 植物の大きさ
  • f : 育てる時の肥料の有無(Cが肥料なし、Tが肥料あり)

です。種子の数について、植物の大きさや肥料との関係を調べます。
個体ごとに飛ばす種子の数の平均値が違うと仮定します。

先にネタバレしてしまうと、種子(y)と植物の大きさ(x)には関係がありますが、種子(y)と肥料の有無(f)は特に関係ありません。
肥料の有無はどちらかというと植物の大きさに関係してるので、確かに直感通りな結果ではあります。

一般化線形モデル

実際に関係を調べてみましょう。
天下りではありますが、fはyと直接関係のないデータなので、yとxの関係だけを調べます。
統計モデルとしてはポアソン回帰を用います。*1

数式で表すと、ある個体iが種子を飛ばす個数は、
{ \displaystyle
p(y_i|\lambda_i) = \frac{\lambda_{i}^{y_i} \exp(-\lambda_i)}{y_i !}
}
に従うとします。ここで{\lambda_i}は種子を飛ばす数の平均値で、

{ \displaystyle
\log(\lambda_i) = \beta_1 + \beta_{2}x_i
}
と表されるとします。この{\beta_1 , \beta_2}について推定します。

コードはこんな感じです。*2

import statsmodels.api as sm
import pandas

# PandasでCSVの読み込み
data = pandas.read_csv("data3a.csv")

# 植物の大きさ
x = sm.add_constant(data.x)
# 飛ばした種子の数
y = data.y

# ポアソン回帰
model = sm.GLM(y, x, family=sm.families.Poisson())
results = model.fit()

# 結果の表示
print(results.summary())

APIの詳細はこちら。
http://statsmodels.sourceforge.net/devel/glm.html
http://statsmodels.sourceforge.net/devel/generated/statsmodels.tools.tools.add_constant.html

一つだけ注意することがあって、sm.add_constant でxに定数列を加えないと、切片({\beta_1})抜きで推定してしまいます。ここだけ何か腑に落ちないです。

結果はこんな感じです。

                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       98
Model Family:                 Poisson   Df Model:                            1
Link Function:                    log   Scale:                             1.0
Method:                          IRLS   Log-Likelihood:                -235.39
Date:                Tue, 13 May 2014   Deviance:                       84.993
Time:                        11:39:50   Pearson chi2:                     83.8
No. Iterations:                     6
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const          1.2917      0.364      3.552      0.000         0.579     2.005
x              0.0757      0.036      2.125      0.034         0.006     0.145
==============================================================================

coefの列が推定値で、const{\beta_1}、xが{\beta_2}に対応しています。

こんな感じで(調べた時間の割には)あっさりとしたコードでポアソン回帰できました。
一番大切な結果の考察などは教科書をご参照あれ。良書です。

参考

私より詳しく解説してらっしゃる方がいました。
https://github.com/Youki/statistical-modeling-for-data-analysis-with-python/blob/master/Chapter03/Chapter03.rst

StatsModels はRのインターフェイスに似せて使えるみたいです。
その際のAPIはこちらです。
Generalized Linear Models (Formula) — statsmodels 0.6.0 documentation

補遺

scikit-learnにももちろん一般化線形モデルが実装されています。
1.1. Generalized Linear Models — scikit-learn 0.14 documentation

紹介した本のサポートページも充実してるので是非。
生態学データ解析 - 本/データ解析のための統計モデリング入門

*1:ポアソン回帰を用いる理由としては、yが離散値だから、平均と分散がだいたい一致してそうだから、などがあります。詳しくは久保先生の本を。

*2:ダウンロードしたデータと同じフォルダにスクリプトファイルがあると仮定しています。