Kubo_book_chap3

3182 days ago by takepwave

Hiroshi TAKEMOTO (take@pwv.co.jp)

Sageでグラフを再現してみよう：データ解析のための統計モデリング入門第３章

この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。

前回に続き、データ解析のための統計モデリング入門 （以下、久保本と書きます） の第３章の例題をSageを使って再現してみます。

数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるの大きな特徴です。 この機会にSageを分析に活用してみてはいかがでしょう。

前準備

最初に必要なライブラリーやパッケージをロードしておきます。新しくなったRUtil.pyも使います。

# Rの必要なライブラリ #r("install.packages('ggplot2')") r('library(ggplot2)') r('library(jsonlite)') # python用のパッケージ import pandas as pd import numpy as np import matplotlib.pyplot as plt from ggplot import * # RUtilにRとPandasのデータフレームを相互に変換する関数を追加 load(DATA + 'RUtil.py') 
       

例題のデータ

３章の例題に使われているデータは、架空の植物の第i番目の個体のサイズ$x_i$と肥料の有無$f_i$が種子数$y_i$ にどう影響するかを調べます。 本の図３．１を以下に引用します。

３章のデータは、ネット公開されており、以下の例でもネットのデータを利用します。

データをロードしたら、summary関数とvar関数を使ってデータの素性を確認します。 ここで、確認することは、以下の3項目です。

  • データの個数と欠損値の有無
  • データの平均
  • データの分散
# 3章のデータをネットから取り込む d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv') 
       

カテゴリ

肥料の有無（施肥$f_i$）は、CとTという文字列で表されており、Rでは因子（factor）と呼んでいますが、 アンケートなどのカテゴリと言った方分かりやすいかも知れません。

PandasのデータフレームもRのData.Frameに劣らず、$y_i$は整数型、$x_i$は実数型、$f_i$は因子型と扱ってくれます。

# どんなデータか調べる d.head() 
        
    y      x  f
0   6   8.31  C
1   6   9.44  C
2   6   9.50  C
3  12   9.07  C
4  10  10.16  C
    y      x  f
0   6   8.31  C
1   6   9.44  C
2   6   9.50  C
3  12   9.07  C
4  10  10.16  C
# 個数と平均、最大、最小を、標準偏差をみる d.describe() 
        
                y           x
count  100.000000  100.000000
mean     7.830000   10.089100
std      2.624881    1.008049
min      2.000000    7.190000
25%      6.000000    9.427500
50%      8.000000   10.155000
75%     10.000000   10.685000
max     15.000000   12.400000
                y           x
count  100.000000  100.000000
mean     7.830000   10.089100
std      2.624881    1.008049
min      2.000000    7.190000
25%      6.000000    9.427500
50%      8.000000   10.155000
75%     10.000000   10.685000
max     15.000000   12.400000

データの可視化

施肥の有無（C, T）別に種子種$y_i$とサイズ$x_i$の関係を散布図で表示してみます。 ggplotを使うと、このようなグラフも簡単に表示することができます。（施肥別の指定は、color='f'で行っています）

# F別の分布をみる ggplot(d, aes(x='x', y='y', color='f')) + geom_point() 
        
<ggplot: (8793139305517)>
<ggplot: (8793139305517)>
GgSave('fig-3.2.png') 
        
Saving 11.0 x 8.0 in image.

Saving 11.0 x 8.0 in image.

私は箱ひげ図には馴染みがなく、最初にみたときには驚きました。 慣れてくるとばらつきを知るには良いプロット方法だと感じるようになりました。

この分布の違いは、後で示す施肥別のヒストグラムをみると分かりやすいと思います。

# 箱ひげプロットはggplotにはないので、Rのggplot2でプロットする # junk = r('d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv")') # dをRに渡す PandaDf2RDf(d, 'd') 
       
graph = preGraph("fig3.3.pdf") r('p <- ggplot(d, aes(x=f, y=y)) + geom_boxplot()') r('plot(p)') postGraph(graph) 
       

可視化から見えてくるもの

施肥別の散布図からは、久保本でも整理してあるとおり、

  • サイズxが増加すると共に種子数yが増えているように思える
  • 施肥fの効果はないように思われる
のように感じます。

このようにデータを整理することが大切です。Sageでは処理と文章が1つのドキュメントにまとめられており、 バージョン管理もされているので、どのようにして結果を導いたのかを後から簡単にフォローすることができます。

施肥の違いをもう少しみてみる

PandasのデータフレームもRのData.Frame同様にデータの絞り込みができます。

施肥のあるもの$f_T$を、以下の様にfのフィールドが'T"のものとして、d_Tにセットすることが1行で記述できます。 

	d_T = d[d['f'] == 'T']
また、Pandasでは、列（Series）をドット(.)で指定することもできるので、以下のようにも記述できます。 
	d_T = d[d.f == 'T']

# 施肥の有無の違いをみる d_T = d[d['f'] == 'T'] print d_T.describe(), d_T.var() 
        
              y          x
count  50.00000  50.000000
mean    7.88000  10.370600
std     2.65453   0.944307
min     3.00000   8.520000
25%     6.00000   9.757500
50%     7.50000  10.225000
75%     9.00000  10.947500
max    15.00000  12.400000 y    7.046531
x    0.891716
dtype: float64
              y          x
count  50.00000  50.000000
mean    7.88000  10.370600
std     2.65453   0.944307
min     3.00000   8.520000
25%     6.00000   9.757500
50%     7.50000  10.225000
75%     9.00000  10.947500
max    15.00000  12.400000 y    7.046531
x    0.891716
dtype: float64

同様に施肥のないもの$f_C$を、抽出します。

平均と分散がほぼ同じであることから、どちらもポアソン分布と仮定しても良さそうです。

d_C = d[d['f'] == 'C'] print d_C.describe(), d_C.var() 
        
               y          x
count  50.000000  50.000000
mean    7.780000   9.807600
std     2.620874   0.999813
min     2.000000   7.190000
25%     6.000000   9.115000
50%     8.000000   9.880000
75%    10.000000  10.412500
max    12.000000  11.930000 y    6.868980
x    0.999627
dtype: float64
               y          x
count  50.000000  50.000000
mean    7.780000   9.807600
std     2.620874   0.999813
min     2.000000   7.190000
25%     6.000000   9.115000
50%     8.000000   9.880000
75%    10.000000  10.412500
max    12.000000  11.930000 y    6.868980
x    0.999627
dtype: float64

列の度数分布は、以下の様にgroupbyでy毎の頻度を求めることができます。

# 施肥なし（C）の度数を調べる d_C.groupby('y').size() 
        
y
2      1
3      2
4      3
5      4
6     10
7      1
8      5
9      8
10     9
11     4
12     3
dtype: int64
y
2      1
3      2
4      3
5      4
6     10
7      1
8      5
9      8
10     9
11     4
12     3
dtype: int64

施肥別のヒストグラム

ggplotを使えば、施肥別のヒストグラムも以下の1行でプロットできます。

facet_wrapでfを指定することで施肥別のヒストグラムを図化できます。

# F別のヒストグラムをみる ggplot(d, aes(x="y")) + geom_histogram(color="grey") + facet_wrap("f") 
        
<ggplot: (8793138564173)>
<ggplot: (8793138564173)>
GgSave('fig-3.0.png') 
        
Saving 11.0 x 8.0 in image.

Saving 11.0 x 8.0 in image.

ポアソン回帰の統計モデル

施肥別の平均、分散からも$個体_i$の種子数$y_i$の確率$p(y_i | \lambda_i)$は、ポアソン分布に従うとして、 以下の式で表します。 $$ p(y_i | \lambda_i) = \frac{\lambda_i^{y_i}exp(-\lambda_i)}{y_i!} $$

$個体_iの平均種子数\lambda_i$の式

久保本の3.4.1では、$個体_iの平均種子数\lambda_i$を以下の様な指数関数で表しています。 $$ \lambda_i = exp( \beta_1 + \beta_2 x_i) $$

ここで、$\beta_1, \beta_2$は、式を決定づけるパラメータです。図3.4では、このパラメータを変えることで、 どうような曲線になるか示しています。これをSageを使って表示すると以下のようになります。

直線よりも指数関数の表現力があるようにみえます。

# fig3.4 p1 = plot(lambda x: exp(-1+0.4*x), [x, -4, 5], rgbcolor='red') p2 = plot(lambda x: exp(-2-0.8*x), [x, -4, 5], rgbcolor='blue') (p1 + p2).show(figsize=5) 
       

線形予測子とリンク関数

glmで使われるリンク関数と線形予測子の関係を整理すると以下の様になります。 線形予測子は、パラメータが線形結合しているもので、ターゲットとなる変数 （ここでは$\lambda_i$）を線形予測子で表したときに、以下の関係がある場合、 $$ \lambda_i = f(線形予測子) $$ リンク関数は、fの逆関数を使って以下の関係にあります。 $$ f^{-1}(\lambda_i) = (線形予測子) $$

ですから指数関数で表された$\lambda_i$のリンク関数は、対数を使って以下の様になります。 $$ log \lambda_i = \beta_1 + \beta_2 x_i $$

ポアソン回帰の実行

久保本では、モデルから推定される予測の良さを重視しており、その指標として対数尤度が最大になることを目安としています。

今回のデータに対する対数尤度(Log-Likelihood)は、以下の様になります。この対数尤度の最大値を追求することになります。 $$ log L(\beta_1, \beta_2) = \sum_i log \frac{\lambda_i^{y_i} exp(\lambda_i)}{y_i !} $$

ここでは、Pythonのstatsmodelsパッケージを使ってGLM（一般化線型モデル）を使ってポアソン回帰を行います。

# statsmodelsを使ってglmを計算します import statsmodels.formula.api as smf import statsmodels.api as sm from scipy.stats.stats import pearsonr 
       

glmの引数

glmの引数について説明すると、

  • fit: 回帰の結果をセットする変数
  • y ~ x: モデル式を線形予測子で表現します
  • data: dでデータフレームを指定します
  • family: sm.families.Poissonで確率分布としてポアソン分布を指定します
  • link: sm.families.links.logで対数を指定します
  • fit(): 最後にfit関数で回帰を計算します

回帰の結果は、summary2関数で表示します。結果は、

  • 対数尤度(Log-Likelihood): -235.39
  • 逸脱度(deviance): 84.993
  • 切片$\beta_1$(Intercept): 1.2917、標準誤差は0.3637
  • xの係数$\beta_2$: 0.0757、標準誤差は0.0356
と求まります。

fit = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit.summary2() 
        
<class 'statsmodels.iolib.summary2.Summary'>
"""
          Results: Generalized linear model
=====================================================
Model:               GLM   AIC:             474.7725 
Link Function:       log   BIC:             -366.3137
Dependent Variable:  y     Log-Likelihood:  -235.39  
No. Observations:    100   Deviance:        84.993   
Df Model:            1     Pearson chi2:    83.8     
Df Residuals:        98    Scale:           1.0000   
Method:              IRLS                            
-----------------------------------------------------
          Coef.  Std.Err.   t    P>|t|  [0.025 0.975]
-----------------------------------------------------
Intercept 1.2917   0.3637 3.5517 0.0004 0.5789 2.0045
x         0.0757   0.0356 2.1251 0.0336 0.0059 0.1454
=====================================================

"""
<class 'statsmodels.iolib.summary2.Summary'>
"""
          Results: Generalized linear model
=====================================================
Model:               GLM   AIC:             474.7725 
Link Function:       log   BIC:             -366.3137
Dependent Variable:  y     Log-Likelihood:  -235.39  
No. Observations:    100   Deviance:        84.993   
Df Model:            1     Pearson chi2:    83.8     
Df Residuals:        98    Scale:           1.0000   
Method:              IRLS                            
-----------------------------------------------------
          Coef.  Std.Err.   t    P>|t|  [0.025 0.975]
-----------------------------------------------------
Intercept 1.2917   0.3637 3.5517 0.0004 0.5789 2.0045
x         0.0757   0.0356 2.1251 0.0336 0.0059 0.1454
=====================================================

"""
# 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットされている print fit.deviance, fit.llf, fit.aic 
        
84.9929964907 -235.38625077 474.77250154
84.9929964907 -235.38625077 474.77250154

回帰結果のプロット

λの予測値は、predict関数で取得することができます（fitのnuにセットされているみたいです）。

施肥別の分布図に予測値を重ね書きします。 予測値をyhatにセットし、ggplotで以下の様にgeom_lineを追加します。

# λの予測値nuをyhatにセット d['yhat'] = fit.predict() 
       
# λの予測値と一緒にプロット ggplot(d, aes(x='x', y='y', color='f')) + geom_point() + geom_line(aes(x='x', y='yhat')) 
        
<ggplot: (8793155860005)>
<ggplot: (8793155860005)>
GgSave('fig-3.7.png') 
        
Saving 11.0 x 8.0 in image.

Saving 11.0 x 8.0 in image.

施肥処理Fをモデルに組み込む

つぎに施肥処理fを使ってポアソン回帰をしてみましょう。xと異なりfには数字では無くT, Cのような文字列（因子型）が入っています。

因子型で線形予測子を作る場合、因子の数より1個少ない数の$d_i$というダミー変数を使って、因子の影響を線形予測子として表します。 $d_iが、f_i=Tの場合1，f_i!=Tの場合0$とします。因子型ではいずれか一つの因子に当てはまる（昔のカーラジオのボタンのように） としているので、$f_iがTでなければCとなるので、因子がTとCの2個の場合ダミー変数は、$d_i$1個で足ります。

もし、肥料が複数の場合には、ダミー変数は、$d_i,A, d_i,B, d_i,C$のように複数作らなくてはなりません。

glmは、とてもお利口なので人がダミー変数を指定する代わりに因子型の変数fを指定するだけで内部的にダミー変数を考慮した計算をしてくれます。

モデル式をy ~ xからy ~ fに変えてポアソン回帰を実行してみます。

回帰の結果は、summary2関数で表示します。結果は、

  • 対数尤度(Log-Likelihood): -237.63
  • 逸脱度(deviance): 89.475
  • 切片$\beta_1$(Intercept): 2.0516、標準誤差は0.0507
  • fの係数$\beta_3$: 0.0128、標準誤差は0.0715
先の例と比べると、対数尤度は-237.63と小さく求まっています。

久保本は、とても丁寧でこの結果をどう考えるかを解説してくれています。 $$ \lambda_i = exp(2.05 + 0.0128) = exp(2.05) exp(0.0128) $$ 施肥を行った場合には、施肥を行わない場合（0.0128の項が０）の1.0128倍に若干増えると予測しています。 $$ \lambda_i = exp(2.05 + 0.0128) = 1.0128 exp(2.05) $$

# 施肥処理Fをモデルに組み込む fit_f= smf.glm('y ~ f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_f.summary2() 
        
<class 'statsmodels.iolib.summary2.Summary'>
"""
           Results: Generalized linear model
=======================================================
Model:                 GLM   AIC:             479.2545 
Link Function:         log   BIC:             -361.8317
Dependent Variable:    y     Log-Likelihood:  -237.63  
No. Observations:      100   Deviance:        89.475   
Df Model:              1     Pearson chi2:    87.1     
Df Residuals:          98    Scale:           1.0000   
Method:                IRLS                            
-------------------------------------------------------
          Coef.  Std.Err.    t    P>|t|   [0.025 0.975]
-------------------------------------------------------
Intercept 2.0516   0.0507 40.4630 0.0000  1.9522 2.1509
f[T.T]    0.0128   0.0715  0.1787 0.8582 -0.1273 0.1529
=======================================================

"""
<class 'statsmodels.iolib.summary2.Summary'>
"""
           Results: Generalized linear model
=======================================================
Model:                 GLM   AIC:             479.2545 
Link Function:         log   BIC:             -361.8317
Dependent Variable:    y     Log-Likelihood:  -237.63  
No. Observations:      100   Deviance:        89.475   
Df Model:              1     Pearson chi2:    87.1     
Df Residuals:          98    Scale:           1.0000   
Method:                IRLS                            
-------------------------------------------------------
          Coef.  Std.Err.    t    P>|t|   [0.025 0.975]
-------------------------------------------------------
Intercept 2.0516   0.0507 40.4630 0.0000  1.9522 2.1509
f[T.T]    0.0128   0.0715  0.1787 0.8582 -0.1273 0.1529
=======================================================

"""

個体のサイズと施肥を合わせたモデル

最後に個体のサイズ$x_i$と施肥$f_i$を合わせたモデルを使ってポアソン回帰を行ってみます。

モデル式は、y ~ x + fとなります。

結果は以下の様になります。 $$ \lambda_i = 1.26 + 0.08 x_i - 0.032 $$

# 説明変数が数量型＋因子型のモデル fit_all = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_all.summary2() 
        
<class 'statsmodels.iolib.summary2.Summary'>
"""
            Results: Generalized linear model
========================================================
Model:                GLM    AIC:              476.5874 
Link Function:        log    BIC:              -361.8936
Dependent Variable:   y      Log-Likelihood:   -235.29  
No. Observations:     100    Deviance:         84.808   
Df Model:             2      Pearson chi2:     83.8     
Df Residuals:         97     Scale:            1.0000   
Method:               IRLS                              
--------------------------------------------------------
           Coef.  Std.Err.    t    P>|t|   [0.025 0.975]
--------------------------------------------------------
Intercept  1.2631   0.3696  3.4172 0.0006  0.5386 1.9876
f[T.T]    -0.0320   0.0744 -0.4302 0.6670 -0.1778 0.1138
x          0.0801   0.0370  2.1620 0.0306  0.0075 0.1527
========================================================

"""
<class 'statsmodels.iolib.summary2.Summary'>
"""
            Results: Generalized linear model
========================================================
Model:                GLM    AIC:              476.5874 
Link Function:        log    BIC:              -361.8936
Dependent Variable:   y      Log-Likelihood:   -235.29  
No. Observations:     100    Deviance:         84.808   
Df Model:             2      Pearson chi2:     83.8     
Df Residuals:         97     Scale:            1.0000   
Method:               IRLS                              
--------------------------------------------------------
           Coef.  Std.Err.    t    P>|t|   [0.025 0.975]
--------------------------------------------------------
Intercept  1.2631   0.3696  3.4172 0.0006  0.5386 1.9876
f[T.T]    -0.0320   0.0744 -0.4302 0.6670 -0.1778 0.1138
x          0.0801   0.0370  2.1620 0.0306  0.0075 0.1527
========================================================

"""

Nullモデルとの比較

ついでなので、λが定数の場合（Nullモデル）のポアソン回帰の結果を以下に示します。

ここで注目するのは、Log-Likelihoodの値-237.64です。 λが定数の場合の対数尤度と自分の考えたモデルの対数尤度を比べることで、 自分のモデルがどの程度よいのか比較できることです。

# Nullモデルの逸脱度とAIC fit_null = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_null.summary2() 
        
<class 'statsmodels.iolib.summary2.Summary'>
"""
           Results: Generalized linear model
======================================================
Model:                GLM   AIC:             477.2864 
Link Function:        log   BIC:             -366.4049
Dependent Variable:   y     Log-Likelihood:  -237.64  
No. Observations:     100   Deviance:        89.507   
Df Model:             0     Pearson chi2:    87.1     
Df Residuals:         99    Scale:           1.0000   
Method:               IRLS                            
------------------------------------------------------
          Coef.  Std.Err.    t    P>|t|  [0.025 0.975]
------------------------------------------------------
Intercept 2.0580   0.0357 57.5862 0.0000 1.9879 2.1280
======================================================

"""
<class 'statsmodels.iolib.summary2.Summary'>
"""
           Results: Generalized linear model
======================================================
Model:                GLM   AIC:             477.2864 
Link Function:        log   BIC:             -366.4049
Dependent Variable:   y     Log-Likelihood:  -237.64  
No. Observations:     100   Deviance:        89.507   
Df Model:             0     Pearson chi2:    87.1     
Df Residuals:         99    Scale:           1.0000   
Method:               IRLS                            
------------------------------------------------------
          Coef.  Std.Err.    t    P>|t|  [0.025 0.975]
------------------------------------------------------
Intercept 2.0580   0.0357 57.5862 0.0000 1.9879 2.1280
======================================================

"""
# 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットされている print fit_null.deviance, fit_null.llf, fit_null.aic 
        
89.5069375696 -237.643221309 477.286442619
89.5069375696 -237.643221309 477.286442619

3章注目のグラフ

「Sageでグラフを再現してみよう」ではいろんなグラフをSageで再現してその本質を理解できるように挑戦しています。

3章の注目のグラフは、図3.9（以下久保本から引用します）の（Ｂ）です。 種子数λがサイズxに対して指数関数で増加するとそのポアソン分布は裾野が広く平らに分布しています。 今回はこのグラフをSageを使って計算してみます。

GLMの結果が示されていないので、$x$が0.5, 1.1, 1.7のyの値をグラフから読み取り、 0.1, 0.5, 3.0としてポアソン分布をプロットしたのが、以下の図です。図３．９（Ｂ）のポアソン分布 とよく似た分布になっています。

# 3章のメインのグラフは、図3.9だと思う。 # 0.5, 1.1, 1.7のyを0.1, 0.5, 3.0としてポアソン分布を表示 # Fig. 2.6 # ポアソン関数p(y, λ)を定義 def p(y, lam): return lam^y*e^(-lam)/factorial(y) g = Graphics() for lam, cls in [(0.1, 'red'), (0.5, 'blue'), (3.0, 'green')]: g = g + list_plot([p(y, lam) for y in range(8)], color=cls) g.show(figsize=5)