#!/usr/bin/env python # coding: utf-8 # # Hiroshi TAKEMOTO # (take.pwave@gmail.com) # #

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

#

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

#

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

#

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

# # # #

前準備

#

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

# # In[1]: get_ipython().run_cell_magic('HTML', '', '\n') # In[2]: # RとPandasのデータフレームを相互に変換する関数を読み込む # Rの必要なライブラリ r('library(ggplot2)') r('library(jsonlite)') # python用のパッケージ import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns get_ipython().run_line_magic('matplotlib', 'inline') # jupyter用のdisplayメソッド from IPython.display import display, Latex, HTML, Math, JSON # sageユーティリティ load('script/sage_util.py') # Rユーティリティ load('script/RUtil.py') # #

例題のデータ

#

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

#

# #

#

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

#

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

#

# # # In[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$は因子型と扱ってくれます。 #

# # In[4]: # どんなデータか調べる d.head() # In[5]: # 個数と平均、最大、最小を、標準偏差をみる d.describe() # #

データの可視化

#

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

# # In[6]: # F別の分布をみる sns.lmplot('x', 'y', data=d, hue='f', fit_reg=False ) plt.show() # #

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

#

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

# # In[7]: # F別の箱ひげ図 sns.boxplot(x='f', y='y', data=d) plt.show() # #

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

#

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

# のように感じます。 #

#

# このようにデータを整理することが大切です。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']	
#
#

# # In[8]: # 施肥の有無の違いをみる d_T = d[d['f'] == 'T'] d_T.describe() # In[9]: d_T.var() # #

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

#

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

# # In[10]: d_C = d[d['f'] == 'C'] d_C.describe() # In[11]: d_C.var() # #

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

# # In[12]: # 施肥なし（C）の度数を調べる d_C.groupby('y').size() # #

施肥別のヒストグラム

#

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

#

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

# # In[13]: # F別のヒストグラムをみる d.hist(column='y', by='f', sharex=True, layout=(2,1), bins=11) plt.show() # #

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

#

# 施肥別の平均、分散からも$個体_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を使って表示すると以下のようになります。 #

#

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

# # In[14]: # 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=4) # #

線形予測子とリンク関数

#

# 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（一般化線型モデル）を使ってポアソン回帰を行います。 #

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

glmの引数

#

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

#

#

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

# と求まります。 #

# # In[16]: fit = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit.summary2() # In[17]: # 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットされている print fit.deviance, fit.llf, fit.aic # #

回帰結果のプロット

#

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

#

# 施肥別の分布図に予測値を重ね書きします。 # 予測値をyhatにセットし、pandasでソートした後、lineplotします。 #

# # In[18]: # λの予測値nuをyhatにセット d['yhat'] = fit.predict() # In[19]: ax = d.sort_values(by=['x']).plot(x='x', y='yhat') d_T.plot(kind='scatter', x='x', y='y', c='r', ax=ax, label='T') d_C.plot(kind='scatter', x='x', y='y', c='b', ax=ax, label='C') plt.legend(loc='upper left') plt.show() # #

施肥処理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関数で表示します。結果は、 #

# 先の例と比べると、対数尤度は-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) # $$ #

# # In[20]: # 施肥処理Fをモデルに組み込む fit_f= smf.glm('y ~ f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_f.summary2() # #

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

#

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

#

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

#

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

# # In[21]: # 説明変数が数量型＋因子型のモデル fit_all = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_all.summary2() # #

Nullモデルとの比較

#

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

#

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

# # In[22]: # Nullモデルの逸脱度とAIC fit_null = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_null.summary2() # In[23]: # 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットされている print fit_null.deviance, fit_null.llf, fit_null.aic # #

3章注目のグラフ

#

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

#

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

#

# #

#

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

# # # In[24]: # 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 float(lam^y*e^(-lam)/factorial(y) ) # データフレームを作成 df = pd.DataFrame( { '0.1': [p(y, 0.1) for y in range(8)], '0.5': [p(y, 0.5) for y in range(8)], '3.0': [p(y, 3.0) for y in range(8)] }) # データフレームのプロット df.plot.bar() plt.show() # In[ ]: