![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2014/04/06からのアクセス回数 5946
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/42/
また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。
前回に続き、<a href="http://www.amazon.co.jp/dp/400006973X/">データ解析のための統計モデリング入門</a> （以下、久保本と書きます） の第5章の例題をSageを使って再現してみます。
私が統計を嫌いになったのは、χ二乗検定が原因です。 どうしてそうなるのかも説明されないまま検定をすることに納得がいかず拒否していたのを覚えています。
久保本はそれをすぱっと解決してくれました。今回の目標はずばり図５．４（以下に久保本から引用）です。
数式処理システム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とでました。
sageへの入力:
# 3章のデータをもう一度読み込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv')
sageへの入力:
# 一定モデル（k=1）でのglm回帰を実行
fit1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()
sageへの入力:
# xモデル
fit2 = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit()
sageへの入力:
# 逸脱度（-2logL） -2*(fit1.llf - fit2.llf)
4.5139410788522696
sageへの入力:
# 切片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
第一種の過誤の検討に専念して、以下の手順で帰無仮説を棄却します。
Sageには、ポアソン分布を生成する関数が無いので、Rのrpois関数を使ってlambda_targetで指定されたλのサンプルを sample_size個生成する関数genSampleを以下の様に定義します。
sageへの入力:
# 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回繰り返し、その分布を求めます。
sageへの入力:
# パラメトリックブートストラップ法
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への入力:
# 逸脱度の差を1000サンプル生成 dd12 = pb(d, 1000)
ブートストラップ法で求まった結果にdescribe（Rのsummaryに相当）を使って分布の概要を把握し、 ggplotのヒストグラムでその分布を表示してみます。
sageへの入力:
dd12Df = pd.DataFrame({'x' : dd12})
dd12Df.describe()
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への入力:
ggplot(dd12Df, aes(x="x")) + geom_histogram(color="grey") + geom_vline(x=4.5)
binwidth defaulted to range/30. Use 'binwidth = x' to adjust this. <ggplot: (13280017)>
sageへの入力:
ggsave('fig-5.4.png', dpi=50)
Saving 11.0 x 8.0 in image.
4.5以上のデータの個数は、37個でp値は、37/1000=0.037となり、95%となるｘの値は 3.90と求まりました。
これで、有意水準5%の検定において帰無仮説（一定モデル）は棄却され、xモデルが採択されました。
sageへの入力:
# 4.5以上のデータの個数は len(dd12Df[dd12Df.x >= 4.5])
37
sageへの入力:
# 95%となるxの値を求める dd12Df.quantile(0.95)
x 3.899619 dtype: float64
残念ながら、statsmodelsにはglmに対するanovaがないため、Rにデータを渡して、計算しました。 p値は、0.034となり、同様に帰無仮説は棄却されました。（Rを使っているので、久保本と同じ結果になるのは当然です。）
sageへの入力:
# χ自乗分布を使かう方法
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への入力:
# 残念ながらstatsmodelsにはglmに対するanovaがないみたい from statsmodels.stats.anova import anova_lm anova_lm(fit1, fit2)
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への入力:
# Rにｄを渡して計算 PandaDf2RDf(d, "d")
sageへの入力:
# 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
皆様のご意見、ご希望をお待ちしております。
