Hiroshi TAKEMOTO (take@pwv.co.jp)

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

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

前回に続き、データ解析のための統計モデリング入門（以下、久保本と書きます）の第5章の例題を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')

尤度比検定がすごい

どんな統計モデルにも使える尤度比検定が、5章のテーマです。ターゲットの統計モデル（xモデル）と一定値を使った意味の無いモデル：帰無仮説の尤度比を使って帰無仮説棄却の可否を判断します。（久保本では「無に帰される」から帰無仮説と説明があり昔の人はよい名前を付けるなぁと感動！）

今回は、3章の例題で一定モデルとｘモデルを比較します。対数尤度は、以下の様になり、 $$ \frac{L^*_1}{L^*_2} = \frac{一定モデルの最大尤度}{ｘモデルの最大尤度} $$ これの対数に-2を掛けた値（逸脱度の差）を使って検定をします。 $$ \Delta D_{1,2} = -2 \times ( log L^*_1 - log L^*_2) $$

逸脱度の差は、以下の計算で4.5とでました。

# 3章のデータをもう一度読み込む d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv')

# 一定モデル（k=1）でのglm回帰を実行 fit1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()

# xモデル fit2 = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()

# 逸脱度（-2logL） -2*(fit1.llf - fit2.llf)

4.5139410788522696

4.5139410788522696

# 切片beta_2の推定値2.058となり、exp(2.058)がd.yの平均とだいたい一致していることを確認 print fit1.summary2() print exp(2.058), d.y.mean()

           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
======================================================

7.83029355218137 7.83

           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
======================================================

7.83029355218137 7.83

帰無仮説を棄却させる手順

第一種の過誤の検討に専念して、以下の手順で帰無仮説を棄却します。

帰無仮説（今回は一定モデル）が正しいものと仮定する
観測データに帰無モデルをあてはめ、$\hat{\beta_1}=2.058$を得て、これが真のモデルとする
上記の真のモデルからサンプルデータを生成し、その度に一定モデルをあてはめ、$\Delta D_{1,2}$を求め、その分布を求める
求まった分布から一定モデル（帰無モデル）とxモデルの逸脱度の差$\Delta D_{1,2} \geq 4.5$となる確率Pを評価する

サンプルの生成

Sageには、ポアソン分布を生成する関数が無いので、Rのrpois関数を使ってlambda_targetで指定されたλのサンプルを sample_size個生成する関数genSampleを以下の様に定義します。

# sageobjを使って変換する方が速いのでgenSample関数を使用する def genSample(lambda_target, sample_size): return np.array(sageobj(r('rpois(%d, %f)'%(sample_size, lambda_target))))

パラメトリックブートストラップ法

パラメトリックブートストラップ法をSageで試してみます。4章のバイアスの計算と異なり、一定モデルを真のモデルとしてサンプルデータを生成し、一定モデルとxモデルにGLMを使って回帰分析を行い、対数尤度llfから逸脱度Dの差を求めます。これをn_bootstrap回繰り返し、その分布を求めます。

# パラメトリックブートストラップ法 def pb(d, n_bootstrap): N = len(d) y_mean = d.y.mean() def sampling(d, y_mean, N): # 帰無仮説を真のモデルとしてサンプルデータを生成 d['rnd'] = genSample(y_mean, N) fit1 = smf.glm('rnd ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit2 = smf.glm('rnd ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() return -2*(fit1.llf - fit2.llf) return [ sampling(d, y_mean, N) for i in range(n_bootstrap)]

# 逸脱度の差を1000サンプル生成 dd12 = pb(d, 1000)

サンプルの結果

ブートストラップ法で求まった結果にdescribe（Rのsummaryに相当）を使って分布の概要を把握し、 ggplotのヒストグラムでその分布を表示してみます。

dd12Df = pd.DataFrame({'x' : dd12}) dd12Df.describe()

                  x
count  1.000000e+03
mean   9.393243e-01
std    1.309587e+00
min    6.144961e-08
25%    1.010703e-01
50%    3.906226e-01
75%    1.314171e+00
max    1.014741e+01

                  x
count  1.000000e+03
mean   9.393243e-01
std    1.309587e+00
min    6.144961e-08
25%    1.010703e-01
50%    3.906226e-01
75%    1.314171e+00
max    1.014741e+01

ggplot(dd12Df, aes(x="x")) + geom_histogram(color="grey") + geom_vline(x=4.5)

<ggplot: (8780298331637)>

<ggplot: (8780298331637)>

GgSave('fig-5.4.png')

Saving 11.0 x 8.0 in image.

Saving 11.0 x 8.0 in image.

仮説の検証

4.5以上のデータの個数は、37個でp値は、37/1000=0.037となり、95%となるｘの値は 3.90と求まりました。

これで、有意水準5%の検定において帰無仮説（一定モデル）は棄却され、xモデルが採択されました。

# 4.5以上のデータの個数は len(dd12Df[dd12Df.x >= 4.5])

# 95%となるxの値を求める dd12Df.quantile(0.95)

x    3.53806
dtype: float64

x    3.53806
dtype: float64

χ自乗分布を使かう方法

残念ながら、statsmodelsにはglmに対するanovaがないため、Rにデータを渡して、計算しました。 p値は、0.034となり、同様に帰無仮説は棄却されました。（Rを使っているので、久保本と同じ結果になるのは当然です。）

# χ自乗分布を使かう方法 fit1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit2 = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()

# 残念ながらstatsmodelsにはglmに対するanovaがないみたい from statsmodels.stats.anova import anova_lm anova_lm(fit1, fit2)

Traceback (click to the left of this block for traceback)
...
AttributeError: 'GLMResults' object has no attribute 'ssr'

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_20.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("IyDmrovlv7XjgarjgYzjgolzdGF0c21vZGVsc+OBq+OBr2dsbeOBq+WvvuOBmeOCi2Fub3Zh44GM44Gq44GE44G/44Gf44GECmZyb20gc3RhdHNtb2RlbHMuc3RhdHMuYW5vdmEgaW1wb3J0IGFub3ZhX2xtCmFub3ZhX2xtKGZpdDEsIGZpdDIp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpBPi3nD/___code___.py", line 4, in <module>
    exec compile(u'anova_lm(fit1, fit2)
  File "", line 1, in <module>
    
  File "/usr/local/sage-6.7/local/lib/python2.7/site-packages/statsmodels-0.6.0-py2.7-linux-x86_64.egg/statsmodels/stats/anova.py", line 353, in anova_lm
    table["ssr"] = map(getattr, args, ["ssr"]*n_models)
  File "/usr/local/sage-6.7/local/lib/python2.7/site-packages/statsmodels-0.6.0-py2.7-linux-x86_64.egg/statsmodels/base/wrapper.py", line 35, in __getattribute__
    obj = getattr(results, attr)
AttributeError: 'GLMResults' object has no attribute 'ssr'

# Rにｄを渡して計算 PandaDf2RDf(d, "d")

# ANOVA(ANalysis Of VAriance)を使って逸脱度を調べる r('fit1 <- glm(y ~ 1, data=d, family=poisson)') r('fit2 <- glm(y ~ x, data=d, family=poisson)') r('anova(fit1, fit2, test = "Chisq")')

Analysis of Deviance Table

Model 1: y ~ 1
Model 2: y ~ x
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
1        99     89.507                       
2        98     84.993  1   4.5139  0.03362 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Analysis of Deviance Table

Model 1: y ~ 1
Model 2: y ~ x
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
1        99     89.507                       
2        98     84.993  1   4.5139  0.03362 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

kubo_book_chap5

3179 days ago by takepwave