- 追加された行はこの色です。
- 削除された行はこの色です。
[[FrontPage]]
#contents
2014/04/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/42/
また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。
* Sageでグラフを再現してみよう：データ解析のための統計モデリング入門第５章 [#j062beaf]
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。
前回に続き、<a href="http://www.amazon.co.jp/dp/400006973X/">データ解析のための統計モデリング入門</a>
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon.co.jp/dp/400006973X/]]
（以下、久保本と書きます）
の第5章の例題をSageを使って再現してみます。
私が統計を嫌いになったのは、χ二乗検定が原因です。
どうしてそうなるのかも説明されないまま検定をすることに納得がいかず拒否していたのを覚えています。
久保本はそれをすぱっと解決してくれました。今回の目標はずばり図５．４（以下に久保本から引用）です。
&ref(Fig_5.4.png);
数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるのが大きな特徴です。
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#e6f6abc9]
最初に必要なライブラリーやパッケージをロードしておきます。
sageへの入力:
#pre{{
# 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')
}}
** 尤度比検定がすごい [#xe547c53]
どんな統計モデルにも使える尤度比検定が、5章のテーマです。ターゲットの統計モデル（xモデル）と一定値を使った意味の無いモデル
：帰無仮説の尤度比を使って帰無仮説棄却の可否を判断します。
（久保本では「無に帰される」から帰無仮説と説明があり昔の人はよい名前を付けるなぁと感動！）
((久保本では「無に帰される」から帰無仮説と説明があり昔の人はよい名前を付けるなぁと感動！))
今回は、3章の例題で一定モデルとｘモデルを比較します。対数尤度は、以下の様になり、
$$
\frac{L^*_1}{L^*_2} = \frac{一定モデルの最大尤度}{ｘモデルの最大尤度}
$$
これの対数に-2を掛けた値（逸脱度の差）を使って検定をします。
$$
\Delta D_{1,2} = -2 \times ( log L^*_1 - log L^*_2)
$$
逸脱度の差は、以下の計算で4.5とでました。
sageへの入力:
#pre{{
# 3章のデータをもう一度読み込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv')
}}
sageへの入力:
#pre{{
# 一定モデル（k=1）でのglm回帰を実行
fit1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()
}}
sageへの入力:
#pre{{
# xモデル
fit2 = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()
}}
sageへの入力:
#pre{{
# 逸脱度（-2logL）
-2*(fit1.llf - fit2.llf)
}}
#pre{{
4.5139410788522696
}}
sageへの入力:
#pre{{
# 切片beta_2の推定値2.058となり、exp(2.058)がd.yの平均とだいたい一致していることを確認
print fit1.summary2()
print exp(2.058), d.y.mean()
}}
#pre{{
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
}}
** 帰無仮説を棄却させる手順 [#mebe844d]
第一種の過誤の検討に専念して、以下の手順で帰無仮説を棄却します。
+ 帰無仮説（今回は一定モデル）が正しいものと仮定する
+ 観測データに帰無モデルをあてはめ、\(\hat{\beta_1}=2.058\)を得て、これが真のモデルとする
+ 上記の真のモデルからサンプルデータを生成し、その度に一定モデルをあてはめ、\(\Delta D_{1,2}\)を求め、その分布を求める
+ 求まった分布から一定モデル（帰無モデル）とxモデルの逸脱度の差\(\Delta D_{1,2} \geq 4.5\)となる確率Pを評価する
*** サンプルの生成 [#s77e1737]
Sageには、ポアソン分布を生成する関数が無いので、Rのrpois関数を使ってlambda_targetで指定されたλのサンプルを
sample_size個生成する関数genSampleを以下の様に定義します。
sageへの入力:
#pre{{
# sageobjを使って変換する方が速いのでgenSample関数を使用する
def genSample(lambda_target, sample_size):
return np.array(sageobj(r('rpois(%d, %f)'%(sample_size, lambda_target))))
}}
*** パラメトリックブートストラップ法 [#w46fc378]
パラメトリックブートストラップ法をSageで試してみます。4章のバイアスの計算と異なり、
一定モデルを真のモデルとしてサンプルデータを生成し、一定モデルとxモデルにGLMを使って
回帰分析を行い、対数尤度llfから逸脱度Dの差を求めます。
これをn_bootstrap回繰り返し、その分布を求めます。
sageへの入力:
#pre{{
# パラメトリックブートストラップ法
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)]
}}
sageへの入力:
#pre{{
# 逸脱度の差を1000サンプル生成
dd12 = pb(d, 1000)
}}
*** サンプルの結果 [#b0f88167]
ブートストラップ法で求まった結果にdescribe（Rのsummaryに相当）を使って分布の概要を把握し、
ggplotのヒストグラムでその分布を表示してみます。
sageへの入力:
#pre{{
dd12Df = pd.DataFrame({'x' : dd12})
dd12Df.describe()
}}
#pre{{
x
count 1.000000e+03
mean 1.016127e+00
std 1.438180e+00
min 2.531954e-08
25% 9.873545e-02
50% 4.598125e-01
75% 1.367469e+00
max 1.048328e+01
[8 rows x 1 columns]
}}
sageへの入力:
#pre{{
ggplot(dd12Df, aes(x="x")) + geom_histogram(color="grey") + geom_vline(x=4.5)
}}
#pre{{
binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
<ggplot: (13280017)>
}}
sageへの入力:
#pre{{
ggsave('fig-5.4.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-5.4.png);
*** 仮説の検証 [#y8954cad]
4.5以上のデータの個数は、37個でp値は、37/1000=0.037となり、95%となるｘの値は
3.90と求まりました。
これで、有意水準5%の検定において帰無仮説（一定モデル）は棄却され、xモデルが採択されました。
sageへの入力:
#pre{{
# 4.5以上のデータの個数は
len(dd12Df[dd12Df.x >= 4.5])
}}
#pre{{
37
}}
sageへの入力:
#pre{{
# 95%となるxの値を求める
dd12Df.quantile(0.95)
}}
#pre{{
x 3.899619
dtype: float64
}}
** χ自乗分布を使かう方法 [#r56bd0ba]
残念ながら、statsmodelsにはglmに対するanovaがないため、Rにデータを渡して、計算しました。
p値は、0.034となり、同様に帰無仮説は棄却されました。（Rを使っているので、久保本と同じ結果になるのは当然です。）
sageへの入力:
#pre{{
# χ自乗分布を使かう方法
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()
}}
#pre{{
}}
sageへの入力:
#pre{{
# 残念ながらstatsmodelsにはglmに対するanovaがないみたい
from statsmodels.stats.anova import anova_lm
anova_lm(fit1, fit2)
}}
#pre{{
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_80.py", line 10, in <module>
exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" +
途中省略
35, in __getattribute__
obj = getattr(results, attr)
AttributeError: 'GLMResults' object has no attribute 'ssr'
}}
sageへの入力:
#pre{{
# Rにｄを渡して計算
PandaDf2RDf(d, "d")
}}
sageへの入力:
#pre{{
# 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")')
}}
#pre{{
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
}}
** コメント [#n4fa83ed]
#vote(おもしろかった,そうでもない,わかりずらい)
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
