kubo_book_chap6

1528 days ago by takepwave

Hiroshi TAKEMOTO (take@pwv.co.jp)

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

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

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

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

前準備

最初に必要なライブラリーやパッケージをロードしておきます。

# 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 * # statsmodelsを使ってglmを計算します import statsmodels.formula.api as smf import statsmodels.api as sm from scipy.stats.stats import pearsonr # RUtilにRとPandasのデータフレームを相互に変換する関数を追加 load(DATA + 'RUtil.py') 
       

上限のあるカウントデータの回帰分析

6章の最初に二項分布を使った回帰分析の例がでてきます。ポアソン分布では上限のないカウントデータでしたが、 上限のある場合には、二項分布を使うのだと説明がありました。 (自然界ではこのように上限のあるカウントデータが多いので、とても参考になりました)

二項分布用のデータについて

二項分布用のデータは、サポートページにあるdata4a.csvです。 これを読み込んで、データの性質と分布を表示します。

直感的に施肥を施す(f=T)ことでサイズxが大きくなっているのが分かります。 また、サイズと種子数の右上がりの関係が見られます。(種子数の上限は8です。)

# 6章のデータを読み込む d = pd.read_csv(DATA+'data4a.csv') d.head() 
       
   N  y      x  f
0  8  1   9.76  C
1  8  6  10.48  C
2  8  5  10.83  C
3  8  6  10.94  C
4  8  1   9.37  C
   N  y      x  f
0  8  1   9.76  C
1  8  6  10.48  C
2  8  5  10.83  C
3  8  6  10.94  C
4  8  1   9.37  C
d.describe() 
       
         N           y           x
count  100  100.000000  100.000000
mean     8    5.080000    9.967200
std      0    2.743882    1.088954
min      8    0.000000    7.660000
25%      8    3.000000    9.337500
50%      8    6.000000    9.965000
75%      8    8.000000   10.770000
max      8    8.000000   12.440000
         N           y           x
count  100  100.000000  100.000000
mean     8    5.080000    9.967200
std      0    2.743882    1.088954
min      8    0.000000    7.660000
25%      8    3.000000    9.337500
50%      8    6.000000    9.965000
75%      8    8.000000   10.770000
max      8    8.000000   12.440000
ggplot(d, aes(x='x', y='y', color='f')) + geom_point() 
       
<ggplot: (8787834169105)>
<ggplot: (8787834169105)>
GgSave('fig-6.2.png') 
       
Saving 11.0 x 8.0 in image.
Saving 11.0 x 8.0 in image.

二項分布

ここで、Sageのプロット関数を使って二項分布の確率分布がどのような形になっているか 久保本の図6.3と同じ条件でブロットしてみます。

二項分布を $$ p(y | N, q) = \binom{N}{y} q^y (1 -q)^{N-y} $$ 以下の様に_p関数で定義します。

Sageでプロットしてみます。Sageではグラフ(Graphics)にグラフを足すことで、重ね合わせのプロットができます。 とても直感的にグラフを書くことができます。(ggplotも同じです)

# 二項分布を定義 def _p(q, y, N): return binomial(N, y)*q^y*(1-q)^(N-y) 
       
g = Graphics() g += list_plot([_p(0.8, y, 8) for y in (0..8)], plotjoined=True, rgbcolor="red", legend_label ="q=0.8") g += list_plot([_p(0.3, y, 8) for y in (0..8)], plotjoined=True, rgbcolor="green", legend_label ="q=0.3") g += list_plot([_p(0.1, y, 8) for y in (0..8)], plotjoined=True, rgbcolor="blue", legend_label ="q=0.1") g.show(figsize=5) 
       

ロジスティック関数

ロジスティック関数 $$ q_i = logistic(z_i) = \frac{1}{1 + exp(-z_i)} $$ の関数を_logisticとして定義し、その分布をSageを使ってプロットします。

ロジスティック関数のような曲線を持つデータの場合には、二項分布の確率分布があると予測して回帰分析を行います。 リンク関数は、ロジスティック関数の逆関数であるロジットリンク関数を指定します。

ロジット関数は、以下の様に表されます(式a)。 $$ logit(q_i) = log \frac{q_i}{1 - q_i} $$

# ロジスティック関数の定義 def _logistic(z): return 1/(1 + exp(-z)) 
       
# ロジスティック曲線 plot(_logistic(x), [x, -6, 6], figsize=5, legend_label='$q=\\frac{1}{1+exp(-z)}$') 
       

パラメータ推定

残念ながら、statmodelsではcbindに相当する二項分布の解析方法が分からず、 Rを使って計算することにしました。

Sageの中からRの関数を使うことができるので、計算を中断することなく進めることができます。

回帰の結果、$\beta_1 = -19.536, \beta_2 = 1.95, \beta_3=2.02$と求まりました。

# statmodelsでcbind(y, N -y)の部分を表現する方法が見つからなかった # 0, 1の場合には、statmodelsでも解析可能です。 # 残念ですが、Rを使って処理します。 PandaDf2RDf(d, "d") 
       
r('glm(cbind(y, N - y) ~ x + f, data=d, family=binomial)') 
       
Call:  glm(formula = cbind(y, N - y) ~ x + f, family = binomial, data =
d)

Coefficients:
(Intercept)            x           fT  
    -19.536        1.952        2.022  

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:	    499.2 
Residual Deviance: 123 	AIC: 272.2
Call:  glm(formula = cbind(y, N - y) ~ x + f, family = binomial, data = d)

Coefficients:
(Intercept)            x           fT  
    -19.536        1.952        2.022  

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:	    499.2 
Residual Deviance: 123 	AIC: 272.2

結果の図化

ggplot2のstat_smoothで結果の曲線を表示しようとしたのですが、上手くできませんでした。

# glmの結果を使って曲線を表示できませんでした!stat_smoothでご容赦 graph = preGraph("fig-6.7.pdf") #r('p <- ggplot(d, aes(x=x, y=y)) + geom_point() + facet_wrap("f") + stat_smooth(method="glm", family="binomial", formula = cbind(y, N - y) ~ x + f)') r('p <- ggplot(d, aes(x=x, y=y)) + geom_point() + facet_wrap("f") + stat_smooth()') r('plot(p)') postGraph(graph) 
       

そこで、Sageのプロット機能を使って結果を表示してみました。

# 代わりにSageでプロットしてみました。 dT = d[d.f == 'T'] dC = d[d.f == 'C'] dT_pts = list_plot(zip(dT.x, dT.y)) dC_pts = list_plot(zip(dC.x, dC.y), rgbcolor='red') dT_plt = plot(_logistic(-19.536+1.952*x + 2.022)*8, [x, 7, 13]) dC_plt = plot(_logistic(-19.536+1.952*x)*8, [x, 7, 13], rgbcolor='red') (dT_pts+dT_plt+dC_pts+dC_plt).show(figsize=5) 
       

結果の解釈

久保本のすごいところは、単に回帰をすることで終わらず、その解釈を説明しているところです。 リンク関数がロジット関数である(式a)であることから $$ log \frac{q_i}{1 - q_i} = \beta_1 + \beta_2 x_i + \beta_3 f_i $$ となり、両辺の指数を取ると、 $$ \frac{q_i}{1 - q_i} = exp(\beta_1 + \beta_2 x_i + \beta_3 f_i) = exp(\beta_1) exp(\beta_2 x_i) exp(\beta_3 f_i) $$ となります。$\frac{q_i}{1 - q_i}$がオッズと呼ばれる量です。

施肥の有無による違いは、exp(2.02)=7.5倍であると推定されました。

割り算を使わない方法はすごい!

久保本の6章で感動したのは、割り算を使わない方法です。 私のようなものは、割り算を使って正規化してしまいますが、 このように割り算をしないで、回帰分析する方法としてオフセット項の使い方を 紹介しています。

オフセット用のデータには、data4b.csvを使います。 このデータは、

  • 調査地iの面積を$A_i$
  • 調査地iの明るさを$x_i$
  • 調査地iにおける植物個体数$y_i$
の記録です。横軸にA, 縦軸にyを取りデータをプロットすると以下の様になります。 面積が大きいほど植物個体数が多く、xが大きい(明るい)程植物個体数が多いことが 見て取れます。

オフセット項は、offset=np.log(オフセットのカラム)のように指定します。ここでは、オフセット値Aのlogを 取った値が使われていますが、これは以下のような理由からです。

人口密度を以下の様に表し、 $$ \frac{平均個体数\lambda_i}{A_i} = 人口密度 $$ 平均個体数と明るさに指数関数的な関係があると仮定すると、 $$ \lambda_i = A_i \times 人口密度 = A_i exp(\beta_1 + \beta_2 x_i) $$ ここで、面積を掛けている部分をlogを使うことで $$ \lambda_i = exp(\beta_1 + \beta_2 x_i + log A_i) $$ と変形することができます。つまり、$log A_i$だけずれた値(オフセット項)、久保本では「げた」をはかせたと考えると よいとありました。

# オフセット用のデータを読み込む d2 = pd.read_csv(DATA+'data4b.csv') d2.head() 
       
    y     x     A
0  57  0.68  10.3
1  64  0.27  15.6
2  49  0.46  10.0
3  64  0.45  14.9
4  82  0.74  14.0
    y     x     A
0  57  0.68  10.3
1  64  0.27  15.6
2  49  0.46  10.0
3  64  0.45  14.9
4  82  0.74  14.0
ggplot(d2, aes(x='A', y='y', color='x')) + geom_point() 
       
<ggplot: (8787834007825)>
<ggplot: (8787834007825)>
GgSave('fig-6.10-A.png') 
       
Saving 11.0 x 8.0 in image.
Saving 11.0 x 8.0 in image.

推定値

推定された値をプロットすると以下のようになります。 データよりもはっきりと明るさによる影響がきれいに分かります。(そのようにモデルを作ったので当たり前なのですが)

fit = smf.glm('y ~ x', data=d2, offset=np.log(d2.A), family=sm.families.Poisson(link=sm.families.links.log)).fit() fit.summary2() d2['pred'] = fit.predict() 
       
# 予測値のプロット ggplot(d2, aes(x='A', y='pred', color='x')) + geom_point() 
       
<ggplot: (8787833253281)>
<ggplot: (8787833253281)>
GgSave('fig-6.10-C.png') 
       
Saving 11.0 x 8.0 in image.
Saving 11.0 x 8.0 in image.

オフセット項を使わないと

オフセット項を使わないと面積による影響が回帰分析に含まれないため、予測された結果は、ぼんやりとしたものになります。

# 人口密度Aをオフセットを使わない場合 fit_2 = smf.glm('y ~ x', data=d2, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_2.summary2() d2['pred2'] = fit_2.predict() 
       
ggplot(d2, aes(x='A', y='pred2', color='x')) + geom_point() 
       
<ggplot: (8787832643061)>
<ggplot: (8787832643061)>
GgSave('fig-6.10-D.png') 
       
Saving 11.0 x 8.0 in image.
Saving 11.0 x 8.0 in image.

正規分布と尤度

Sageにも確率分布を扱う関数RealDistirubtionがVer.5から導入されました。(まだ一部の関数のサポートですが) これを使って、正規分布の確率密度関数をプロットしていました。

# Sage ver.5から導入されたRealDistributionを使う T1 = RealDistribution('gaussian', 1, seed=100) T2 = RealDistribution('gaussian', 3, seed=100) 
       
T1_plt = plot(lambda x : T1.distribution_function(x), [x, -4, 4], rgbcolor="red") T2_plt = plot(lambda x : T2.distribution_function(x), [x, -4, 4], rgbcolor="green") T3_plt = plot(lambda x : T1.distribution_function(x-2), [x, -4, 4], rgbcolor="blue") (T1_plt + T2_plt + T3_plt).show(figsize=5) 
       

確率密度関数から確率を得る

確率密度関数から指定された範囲に含まれる確率は、その範囲の面積になります。(その範囲で積分した値) cum_distribution_functionは、マイナス無限大から指定された値までの確立を計算するので、 この関数の差を計算することで、指定範囲の確率(面積)を求めることができます。

# xが1.2から1.8となる確率は、cum_distribution_functionの差で計算で求まる T1.cum_distribution_function(1.8) - T1.cum_distribution_function(1.2) 
       
0.07913935110878245
0.07913935110878245

正規分布の尤度

久保本に従って正規分布の尤度をと最尤尤度を求めてみます。 $y_iがy_i - 0.5\Delta y \leq y \leq y_i + 0.5\Delta y$である確率は、 $$ \begin{eqnarray} L(\mu, \sigma) & = & \prod_i p(y_i | \mu, \sigma) \Delta y \\ & = & \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} exp\left\{ \frac{(y_i - \mu)^2}{2\sigma^2} \right\} \Delta y \\ & = & \prod_i (2 \pi \sigma^2)^{-\frac{1}{2}} exp\left\{ \frac{(y_i - \mu)^2}{2\sigma^2} \right\} \Delta y \end{eqnarray} $$ 両辺をlogを取って $$ \begin{eqnarray} log L(\mu, \sigma) & = & -\frac{1}{2} \sum_i log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 + \sum_i log(\Delta y) \\ & = & -\frac{1}{2} N log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 + N log(\Delta y) \end{eqnarray} $$ $\Delta y$の項を無視すると、 $$ log L(\mu, \sigma) = -\frac{1}{2} N log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 $$ となります。これを$\mu$で偏微分し、最小となるのは、$y = \mu$と求まります。(最小自乗法の解と同じ)

ガンマ分布

どんなときにガンマ分布を使うのかGoogleでみると、正の値で、寿命や待ち時間など経過時間に無関係一定値の場合に、 使われるとありました。

ガンマ分布をSageでプロットしてみます。

# Rのdgammaを使用 def _gamma(y, shape, rate): return sageobj(r('dgamma(%s, %s, %s)' %(y, shape, rate))) 
       
g = Graphics() g += plot(lambda x: _gamma(x, 1, 1), [x, 0, 4], rgbcolor="red", legend_label ="r=s=1") g += plot(lambda x: _gamma(x, 5, 5), [x, 0, 4], rgbcolor="green", legend_label ="r=s=5") #g += plot(lambda x: _gamma(x, 0.5, 0.5), [x, 0, 4], rgbcolor="blue", legend_label ="r=s=0.1") g.show(figsize=5) 
       

ガンマ分布を使った回帰分析

例題では、ある個体の花の重量を$y_i$が平均$\mu_i$のガンマ分布にしたがっていると仮定します。 $$ \mu_i = A x_i^b $$ で、A=exp(a)とすると、 $$ \mu_i = exp(a)x_i^b = exp(a + log x_i^b) = exp(a + b log x_i) $$ と表され、リンク関数は、expの逆関数のlogとなり、 $$ log \mu_i = a + b log x_i $$ となります。

ガンマ分布の例題データは、RData形式なので、Rで解析します。family=Gamma(link="log")でリンク関数logのガンマ分布を指定します。

glmを使うと複雑な分布の回帰が簡単に計算できます。しかし、どのモデルが良いかは、 プロット結果を見ながら判断するのがよいと実感しました。人間の目で分からない違いは、尤度やAICを使うのが現実できなのでは?

# ガンマ分布用のデータを使って回帰分析 r('load("%s")'%(DATA + 'd.RData')) r('glm(y ~ log(x), family=Gamma(link="log"), data=d)') 
       
Call:  glm(formula = y ~ log(x), family = Gamma(link = "log"), data = d)

Coefficients:
(Intercept)       log(x)  
    -1.0403       0.6833  

Degrees of Freedom: 49 Total (i.e. Null);  48 Residual
Null Deviance:	    35.37 
Residual Deviance: 17.25 	AIC: -110.9
Call:  glm(formula = y ~ log(x), family = Gamma(link = "log"), data = d)

Coefficients:
(Intercept)       log(x)  
    -1.0403       0.6833  

Degrees of Freedom: 49 Total (i.e. Null);  48 Residual
Null Deviance:	    35.37 
Residual Deviance: 17.25 	AIC: -110.9
# 結果の図化 graph = preGraph("fig-6.13.pdf") r('p <- ggplot(data=d, aes(x=x, y=y)) + geom_point() + stat_smooth(method=glm, formula=y ~ log(x), family=Gamma(link="log"))') r('plot(p)') postGraph(graph)