sage/データ解析のための統計モデリング入門勉強ノート第7章
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2014/04/27からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/44/
また、私の公開しているSageのサーバ(http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう:データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
(以下、久保本と書きます)
の第7章の例題をSageを使って再現してみます。
久保本もようやくテーマの個体差が大きい場合の例題に入りま...
&ref(Fig7.10.png);
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#ke51dbc9]
最初に必要なライブラリーやパッケージをロードしておきます。
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')
}}
** 個体差の大きなデータ [#c7e8434b]
7章の例題は、個体\(x_i\)で観測された以下のデータを含みま...
- 調査種子数\(N_i = 8\)
- 生存種子数\(y_i\)
- 葉数\(x_i\)
- 個体識別番号id
*** データの特徴と可視化 [#lee96c89]
いつもの通りデータを読み込み、その内容、統計情報をチェッ...
sageへの入力:
#pre{{
# 7章のデータを読み込む
d = pd.read_csv(DATA+'data.csv')
d.head()
}}
#pre{{
N y x id
0 8 0 2 1
1 8 1 2 2
2 8 2 2 3
3 8 4 2 4
4 8 1 2 5
[5 rows x 4 columns]
}}
sageへの入力:
#pre{{
d.describe()
}}
#pre{{
N y x id
count 100 100.000000 100.000000 100.000000
mean 8 3.810000 4.000000 50.500000
std 0 3.070534 1.421338 29.011492
min 8 0.000000 2.000000 1.000000
25% 8 1.000000 3.000000 25.750000
50% 8 3.000000 4.000000 50.500000
75% 8 7.000000 5.000000 75.250000
max 8 8.000000 6.000000 100.000000
[8 rows x 4 columns]
}}
*** 同じ点に重なるデータの可視化 [#pd150b89]
プロットデータを見ると、とても100個のデータがあるよう...
これは、種子数・葉数が整数で同じ点に複数のデータが重なっ...
sageへの入力:
#pre{{
# python版のggplotでは、jitter
ggplot(d, aes(x='x', y='y')) + geom_point()
}}
#pre{{
<ggplot: (39068529)>
}}
sageへの入力:
#pre{{
ggsave('fig-7.2.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-7.2.png);
残念ながら、python版のggplotでは同じ点を重ならないように...
Rのggplot2を使ってプロットしてみました。
これで、データのばらつきが正しく把握できました。
sageへの入力:
#pre{{
# Rにdを渡す
PandaDf2RDf(d, "d")
}}
sageへの入力:
#pre{{
# jitterを使って同じポイントに重ならないようにする
graph = preGraph("fig7.2-R.png")
r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=po...
r('plot(p)')
postGraph(graph)
}}
&ref(fig7.2-R.png);
*** GLMを使った分析 [#c7db83b0]
このデータを6章と同じようにGLMで分析してみます。
予測値は直線となっており、6章のようなロジスティック曲線...
sageへの入力:
#pre{{
# 通常のGLMを使って解析
r('fit <- glm(cbind(y, N - y) ~ x , data=d, family=binomi...
r('summary(fit)')
}}
#pre{{
Call:
glm(formula = cbind(y, N - y) ~ x, family = binomial, dat...
Deviance Residuals:
Min 1Q Median 3Q Max
-4.4736 -2.1182 -0.5505 2.3097 4.0966
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1487 0.2372 -9.057 <2e-16 ***
x 0.5104 0.0556 9.179 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ...
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 607.42 on 99 degrees of freedom
Residual deviance: 513.84 on 98 degrees of freedom
AIC: 649.61
Number of Fisher Scoring iterations: 4
}}
sageへの入力:
#pre{{
graph = preGraph("fig7.3-R.png")
r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=po...
r('plot(p)')
postGraph(graph)
}}
&ref(fig7.3-R.png);
久保本に従い、\(x_i = 4\)のデータを抽出し、そのヒストグラ...
二項分布を表示してみましょう。
以下の様にPandasの抽出機能を使ってx=4のデータを取り出し、...
これにロジスティック_logisticと二項分布_pの2つの関数を定...
sageへの入力:
#pre{{
# x = 4の種子数の分布
tbl4 = d[d.x==4].groupby('y').size()
tbl4_lst = zip(tbl4.index, tbl4.tolist())
print tbl4
}}
#pre{{
y
0 3
1 1
2 4
3 2
4 1
5 1
6 2
7 3
8 3
dtype: int64
}}
sageへの入力:
#pre{{
# ロジスティック関数の定義
def _logistic(z):
return 1/(1 + exp(-z))
}}
sageへの入力:
#pre{{
# 二項分布を定義
def _p(q, y, N):
return binomial(N, y)*q^y*(1-q)^(N-y)
}}
GLMで求まったモデルを使って、x=4でのqの値をロジスティック...
x=4での種子数分布(青点)と推定された二項分布(_p関数の値...
推定された種子数になります。
このグラフからは、推定された種子数は、観測された種子数を...
久保本の説明では、\(x_i = 4\)の生存種子数の平均と分散を求...
分散は8.37と大きく、二項分布の分散\(N p (1-p)\)から期待さ...
比べて4倍ほど大きな値となっています。久保本ではこの状態...
sageへの入力:
#pre{{
# x=4での生存確率qを求める
q = _logistic(-2.15+0.51*4)
print q
}}
#pre{{
0.472527695655406
}}
sageへの入力:
#pre{{
# 図7.3(B)をプロット
lst_plt = list_plot(tbl4_lst)
prd_plt = list_plot([_p(q, y, 8)*tbl4.sum() for y in (0.....
(lst_plt + prd_plt).show(figsize=5)
}}
&ref(sage0.png);
sageへの入力:
#pre{{
# x = 4 の平均と分散を求める(平均に対して分散が大きいこ...
d4 = d[d.x==4]
print d4.y.mean(), d4.y.var()
}}
#pre{{
4.05 8.36578947368
}}
sageへの入力:
#pre{{
# 二項分布から期待されている分散は2なので、8.4はかなり大...
N = 8
vr = N*q*(1-q)
print vr
}}
#pre{{
1.99396217995198
}}
** 一般化線形混合分布(GLMM)を使った分析 [#o3efae86]
個体差を考慮した分析手法として一般化線形混合分布(GLMM)...
ここでは、個体\(i\)の個体差を\(r_i\)とし、以下の様なリン...
$$
logit(q_i) = \beta_1 + \beta_2 x_i + r_i
$$
ここで\(r_i\)をどのような分布にするかですが、久保本では平...
しています。(統計モデルとして便利だからと理由です)
glmmMLを使った分析では、個体毎に異なる独立パラメータをclu...
ここでは、idが個体を識別する指標として使われています。
sageへの入力:
#pre{{
#r("install.packages('glmmML')")
r('library(glmmML)')
}}
#pre{{
[1] "glmmML" "jsonlite" "ggplot2" "stats" "gra...
[9] "methods" "base"
}}
sageへの入力:
#pre{{
# python版ではない、R一般線形混合モデルGLMMを使う
# 個体差は、平均0の分散sの正規分布に従うと仮定して先の問...
r('fit <- glmmML(cbind(y, N - y) ~ x, data = d, family = ...
}}
#pre{{
Call: glmmML(formula = cbind(y, N - y) ~ x, family = bin...
coef se(coef) z Pr(>|z|)
(Intercept) -4.190 0.8777 -4.774 1.81e-06
x 1.005 0.2075 4.843 1.28e-06
Scale parameter in mixing distribution: 2.408 gaussian
Std. Error: 0.2202
LR p-value for H_0: sigma = 0: 2.136e-55
Residual deviance: 269.4 on 97 degrees of freedom AI...
}}
*** 結果の可視化 [#lc80752e]
GLMとの違いは、分析結果から予測される二項分布を見ると明ら...
GLMのようなfitted.valuesがGLMMにはないので、分析結果の二...
stat_function使ってプロットしてみました。
きれいに二項分布が表示されました。(これが図7.10の(A))
sageへの入力:
#pre{{
# 結果からロジスティック関数をmyfunc <- 1/(1 + exp(-(-4.1...
r('myfunc <- function(z) { 8*(1/(1 + exp(-(-4.190 + 1.005...
# 結果をプロット
graph = preGraph("fig7.10-R.png")
r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=po...
r('plot(p)')
postGraph(graph)
}}
&ref(fig7.10-R.png);
*** \(r_i\)の分布 [#j0c4243c]
もう一つGLMMから求められた\(r_i\)の正規分布をプロットする...
-10, 10ではほぼ0になっています。
sageへの入力:
#pre{{
# σ=2.4の正規分布
T = RealDistribution('gaussian', 2.4, seed=100)
plot(lambda x : T.distribution_function(x), [x, -10, 10],...
}}
&ref(sage1.png);
*** GLMMから求められた尤度 [#w8626700]
個体毎の尤度\(L_i\)は、以下の様に\(r_i\)を積分して求める...
$$
L_i = \int_{-\infty}^{\infty} p(y_i | \beta_1, \beta_2, r...
$$
GLMMの結果からqを求める関数_qを定義し、上記の積分をSageの...
作りました。rの積分範囲は、sの分布から-10から10としました。
チェックのために、久保本図7.8にあるx=4, r={-2.20, -0.60, ...
y={0, 4, 7}の尤度\(L_i\)を出力してみました。
最後に、Liに20(個体数)を掛けた分布をx=4のヒストグラムに...
きれいに、図7.10(B)の通り、生存種子数を表現できているのが...
このような積分を含む処理でも数式処理システムSageを使えば...
ご理解頂けたと思います。
sageへの入力:
#pre{{
# 回帰分析の結果からqを求める
def _q(x, r):
return _logistic(-4.190 + 1.005*x + r)
# 指定されたyの尤度を求める(上記分布から積分範囲は-10から...
def _Li(y):
return numerical_integral(lambda r: _p(_q(4, r), y, 8...
}}
sageへの入力:
#pre{{
_q(4, -2.2), _q(4, -0.6), _q(4, 1.0), _q(4, 2.6)
}}
#pre{{
(0.0854891394348065, 0.316479106263684, 0.696354929823834...
}}
sageへの入力:
#pre{{
_Li(0), _Li(4), _Li(7)
}}
#pre{{
(0.18814448426796998, 0.07913801041038665, 0.108408446409...
}}
sageへの入力:
#pre{{
prd2_plt = list_plot([_Li(x)*20 for x in range(9)], plotj...
(lst_plt + prd2_plt).show(figsize=5)
}}
&ref(fig7.10_B.png);
** コメント [#ca478220]
#vote(おもしろかった,そうでもない,わかりずらい)
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2014/04/27からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/44/
また、私の公開しているSageのサーバ(http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう:データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
(以下、久保本と書きます)
の第7章の例題をSageを使って再現してみます。
久保本もようやくテーマの個体差が大きい場合の例題に入りま...
&ref(Fig7.10.png);
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#ke51dbc9]
最初に必要なライブラリーやパッケージをロードしておきます。
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')
}}
** 個体差の大きなデータ [#c7e8434b]
7章の例題は、個体\(x_i\)で観測された以下のデータを含みま...
- 調査種子数\(N_i = 8\)
- 生存種子数\(y_i\)
- 葉数\(x_i\)
- 個体識別番号id
*** データの特徴と可視化 [#lee96c89]
いつもの通りデータを読み込み、その内容、統計情報をチェッ...
sageへの入力:
#pre{{
# 7章のデータを読み込む
d = pd.read_csv(DATA+'data.csv')
d.head()
}}
#pre{{
N y x id
0 8 0 2 1
1 8 1 2 2
2 8 2 2 3
3 8 4 2 4
4 8 1 2 5
[5 rows x 4 columns]
}}
sageへの入力:
#pre{{
d.describe()
}}
#pre{{
N y x id
count 100 100.000000 100.000000 100.000000
mean 8 3.810000 4.000000 50.500000
std 0 3.070534 1.421338 29.011492
min 8 0.000000 2.000000 1.000000
25% 8 1.000000 3.000000 25.750000
50% 8 3.000000 4.000000 50.500000
75% 8 7.000000 5.000000 75.250000
max 8 8.000000 6.000000 100.000000
[8 rows x 4 columns]
}}
*** 同じ点に重なるデータの可視化 [#pd150b89]
プロットデータを見ると、とても100個のデータがあるよう...
これは、種子数・葉数が整数で同じ点に複数のデータが重なっ...
sageへの入力:
#pre{{
# python版のggplotでは、jitter
ggplot(d, aes(x='x', y='y')) + geom_point()
}}
#pre{{
<ggplot: (39068529)>
}}
sageへの入力:
#pre{{
ggsave('fig-7.2.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-7.2.png);
残念ながら、python版のggplotでは同じ点を重ならないように...
Rのggplot2を使ってプロットしてみました。
これで、データのばらつきが正しく把握できました。
sageへの入力:
#pre{{
# Rにdを渡す
PandaDf2RDf(d, "d")
}}
sageへの入力:
#pre{{
# jitterを使って同じポイントに重ならないようにする
graph = preGraph("fig7.2-R.png")
r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=po...
r('plot(p)')
postGraph(graph)
}}
&ref(fig7.2-R.png);
*** GLMを使った分析 [#c7db83b0]
このデータを6章と同じようにGLMで分析してみます。
予測値は直線となっており、6章のようなロジスティック曲線...
sageへの入力:
#pre{{
# 通常のGLMを使って解析
r('fit <- glm(cbind(y, N - y) ~ x , data=d, family=binomi...
r('summary(fit)')
}}
#pre{{
Call:
glm(formula = cbind(y, N - y) ~ x, family = binomial, dat...
Deviance Residuals:
Min 1Q Median 3Q Max
-4.4736 -2.1182 -0.5505 2.3097 4.0966
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1487 0.2372 -9.057 <2e-16 ***
x 0.5104 0.0556 9.179 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ...
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 607.42 on 99 degrees of freedom
Residual deviance: 513.84 on 98 degrees of freedom
AIC: 649.61
Number of Fisher Scoring iterations: 4
}}
sageへの入力:
#pre{{
graph = preGraph("fig7.3-R.png")
r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=po...
r('plot(p)')
postGraph(graph)
}}
&ref(fig7.3-R.png);
久保本に従い、\(x_i = 4\)のデータを抽出し、そのヒストグラ...
二項分布を表示してみましょう。
以下の様にPandasの抽出機能を使ってx=4のデータを取り出し、...
これにロジスティック_logisticと二項分布_pの2つの関数を定...
sageへの入力:
#pre{{
# x = 4の種子数の分布
tbl4 = d[d.x==4].groupby('y').size()
tbl4_lst = zip(tbl4.index, tbl4.tolist())
print tbl4
}}
#pre{{
y
0 3
1 1
2 4
3 2
4 1
5 1
6 2
7 3
8 3
dtype: int64
}}
sageへの入力:
#pre{{
# ロジスティック関数の定義
def _logistic(z):
return 1/(1 + exp(-z))
}}
sageへの入力:
#pre{{
# 二項分布を定義
def _p(q, y, N):
return binomial(N, y)*q^y*(1-q)^(N-y)
}}
GLMで求まったモデルを使って、x=4でのqの値をロジスティック...
x=4での種子数分布(青点)と推定された二項分布(_p関数の値...
推定された種子数になります。
このグラフからは、推定された種子数は、観測された種子数を...
久保本の説明では、\(x_i = 4\)の生存種子数の平均と分散を求...
分散は8.37と大きく、二項分布の分散\(N p (1-p)\)から期待さ...
比べて4倍ほど大きな値となっています。久保本ではこの状態...
sageへの入力:
#pre{{
# x=4での生存確率qを求める
q = _logistic(-2.15+0.51*4)
print q
}}
#pre{{
0.472527695655406
}}
sageへの入力:
#pre{{
# 図7.3(B)をプロット
lst_plt = list_plot(tbl4_lst)
prd_plt = list_plot([_p(q, y, 8)*tbl4.sum() for y in (0.....
(lst_plt + prd_plt).show(figsize=5)
}}
&ref(sage0.png);
sageへの入力:
#pre{{
# x = 4 の平均と分散を求める(平均に対して分散が大きいこ...
d4 = d[d.x==4]
print d4.y.mean(), d4.y.var()
}}
#pre{{
4.05 8.36578947368
}}
sageへの入力:
#pre{{
# 二項分布から期待されている分散は2なので、8.4はかなり大...
N = 8
vr = N*q*(1-q)
print vr
}}
#pre{{
1.99396217995198
}}
** 一般化線形混合分布(GLMM)を使った分析 [#o3efae86]
個体差を考慮した分析手法として一般化線形混合分布(GLMM)...
ここでは、個体\(i\)の個体差を\(r_i\)とし、以下の様なリン...
$$
logit(q_i) = \beta_1 + \beta_2 x_i + r_i
$$
ここで\(r_i\)をどのような分布にするかですが、久保本では平...
しています。(統計モデルとして便利だからと理由です)
glmmMLを使った分析では、個体毎に異なる独立パラメータをclu...
ここでは、idが個体を識別する指標として使われています。
sageへの入力:
#pre{{
#r("install.packages('glmmML')")
r('library(glmmML)')
}}
#pre{{
[1] "glmmML" "jsonlite" "ggplot2" "stats" "gra...
[9] "methods" "base"
}}
sageへの入力:
#pre{{
# python版ではない、R一般線形混合モデルGLMMを使う
# 個体差は、平均0の分散sの正規分布に従うと仮定して先の問...
r('fit <- glmmML(cbind(y, N - y) ~ x, data = d, family = ...
}}
#pre{{
Call: glmmML(formula = cbind(y, N - y) ~ x, family = bin...
coef se(coef) z Pr(>|z|)
(Intercept) -4.190 0.8777 -4.774 1.81e-06
x 1.005 0.2075 4.843 1.28e-06
Scale parameter in mixing distribution: 2.408 gaussian
Std. Error: 0.2202
LR p-value for H_0: sigma = 0: 2.136e-55
Residual deviance: 269.4 on 97 degrees of freedom AI...
}}
*** 結果の可視化 [#lc80752e]
GLMとの違いは、分析結果から予測される二項分布を見ると明ら...
GLMのようなfitted.valuesがGLMMにはないので、分析結果の二...
stat_function使ってプロットしてみました。
きれいに二項分布が表示されました。(これが図7.10の(A))
sageへの入力:
#pre{{
# 結果からロジスティック関数をmyfunc <- 1/(1 + exp(-(-4.1...
r('myfunc <- function(z) { 8*(1/(1 + exp(-(-4.190 + 1.005...
# 結果をプロット
graph = preGraph("fig7.10-R.png")
r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=po...
r('plot(p)')
postGraph(graph)
}}
&ref(fig7.10-R.png);
*** \(r_i\)の分布 [#j0c4243c]
もう一つGLMMから求められた\(r_i\)の正規分布をプロットする...
-10, 10ではほぼ0になっています。
sageへの入力:
#pre{{
# σ=2.4の正規分布
T = RealDistribution('gaussian', 2.4, seed=100)
plot(lambda x : T.distribution_function(x), [x, -10, 10],...
}}
&ref(sage1.png);
*** GLMMから求められた尤度 [#w8626700]
個体毎の尤度\(L_i\)は、以下の様に\(r_i\)を積分して求める...
$$
L_i = \int_{-\infty}^{\infty} p(y_i | \beta_1, \beta_2, r...
$$
GLMMの結果からqを求める関数_qを定義し、上記の積分をSageの...
作りました。rの積分範囲は、sの分布から-10から10としました。
チェックのために、久保本図7.8にあるx=4, r={-2.20, -0.60, ...
y={0, 4, 7}の尤度\(L_i\)を出力してみました。
最後に、Liに20(個体数)を掛けた分布をx=4のヒストグラムに...
きれいに、図7.10(B)の通り、生存種子数を表現できているのが...
このような積分を含む処理でも数式処理システムSageを使えば...
ご理解頂けたと思います。
sageへの入力:
#pre{{
# 回帰分析の結果からqを求める
def _q(x, r):
return _logistic(-4.190 + 1.005*x + r)
# 指定されたyの尤度を求める(上記分布から積分範囲は-10から...
def _Li(y):
return numerical_integral(lambda r: _p(_q(4, r), y, 8...
}}
sageへの入力:
#pre{{
_q(4, -2.2), _q(4, -0.6), _q(4, 1.0), _q(4, 2.6)
}}
#pre{{
(0.0854891394348065, 0.316479106263684, 0.696354929823834...
}}
sageへの入力:
#pre{{
_Li(0), _Li(4), _Li(7)
}}
#pre{{
(0.18814448426796998, 0.07913801041038665, 0.108408446409...
}}
sageへの入力:
#pre{{
prd2_plt = list_plot([_Li(x)*20 for x in range(9)], plotj...
(lst_plt + prd2_plt).show(figsize=5)
}}
&ref(fig7.10_B.png);
** コメント [#ca478220]
#vote(おもしろかった,そうでもない,わかりずらい)
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc