FrontPage
2014/04/27からのアクセス回数
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
また、私の公開しているSageのサーバ( )にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。
前回に続き、 http://www.amazon.co.jp/dp/400006973X/ (以下、久保本と書きます) の第7章の例題をSageを使って再現してみます。
久保本もようやくテーマの個体差が大きい場合の例題に入りました。今回のグラフは図7.10です。以下に久保本から引用します。
数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるのが大きな特徴です。 この機会にSageを分析に活用してみてはいかがでしょう。
最初に必要なライブラリーやパッケージをロードしておきます。
sageへの入力:
# Rの必要なライブラリ
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') }}
7章の例題は、個体\(x_i\)で観測された以下のデータを含みます。
いつもの通りデータを読み込み、その内容、統計情報をチェックします。
sageへの入力:
# 7章のデータを読み込む
d = pd.read_csv(DATA+'data.csv') d.head() }}
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への入力:
d.describe() }}
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] }}
プロットデータを見ると、とても100個のデータがあるようには見えません。 これは、種子数・葉数が整数で同じ点に複数のデータが重なっているからです。
sageへの入力:
# python版のggplotでは、jitter
ggplot(d, aes(x='x', y='y')) + geom_point() }}
ggplot: (39068529)> }}
sageへの入力:
ggsave('fig-7.2.png', dpi=50) }}
Saving 11.0 x 8.0 in image. }}
残念ながら、python版のggplotでは同じ点を重ならないようにするjitter機能が使えないため、 Rのggplot2を使ってプロットしてみました。
これで、データのばらつきが正しく把握できました。
sageへの入力:
# Rにdを渡す
PandaDf2RDf(d, "d") }}
sageへの入力:
# jitterを使って同じポイントに重ならないようにする
graph = preGraph("fig7.2-R.png") r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=position_jitter(width=.2, height=0), shape=21, size=2.5)') r('plot(p)') postGraph(graph) }}
このデータを6章と同じようにGLMで分析してみます。 予測値は直線となっており、6章のようなロジスティック曲線にはなっていません。
sageへの入力:
# 通常のGLMを使って解析
r('fit <- glm(cbind(y, N - y) ~ x , data=d, family=binomial)') r('summary(fit)') }}
Call: glm(formula = cbind(y, N - y) ~ x, family = binomial, data = d)
Deviance Residuals:
Coefficients:
(Intercept) -2.1487 0.2372 -9.057 <2e-16 *** x 0.5104 0.0556 9.179 <2e-16 ***
(Dispersion parameter for binomial family taken to be 1)
Residual deviance: 513.84 on 98 degrees of freedom AIC: 649.61
Number of Fisher Scoring iterations: 4 }}
sageへの入力:
graph = preGraph("fig7.3-R.png") r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=position_jitter(width=.2, height=0), shape=21, size=2.5) + geom_line(aes(x=x, y=8*fit$fitted.values))') r('plot(p)') postGraph(graph) }}
久保本に従い、\(x_i = 4\)のデータを抽出し、そのヒストグラムとGLMから推定された 二項分布を表示してみましょう。
以下の様にPandasの抽出機能を使ってx=4のデータを取り出し、頻度を求めます。 これにロジスティック_logisticと二項分布_pの2つの関数を定義します。
sageへの入力:
# x = 4の種子数の分布
tbl4 = d[d.x==4].groupby('y').size() tbl4_lst = zip(tbl4.index, tbl4.tolist()) print tbl4 }}
y 0 3 1 1 2 4 3 2 4 1 5 1 6 2 7 3 8 3 dtype: int64 }}
sageへの入力:
# ロジスティック関数の定義
def _logistic(z):
}}
sageへの入力:
# 二項分布を定義
def _p(q, y, N):
}}
GLMで求まったモデルを使って、x=4でのqの値をロジスティック関数から求めます。
x=4での種子数分布(青点)と推定された二項分布(_p関数の値)にデータ数(tbl4.sum()=20)を掛けた値が、 推定された種子数になります。
このグラフからは、推定された種子数は、観測された種子数を説明しているように思えません。 久保本の説明では、\(x_i = 4\)の生存種子数の平均と分散を求めてみると、平均4.05に対して、 分散は8.37と大きく、二項分布の分散\(N p (1-p)\)から期待される値\(8\times0.5\times(1-0.5)=2\)と 比べて4倍ほど大きな値となっています。久保本ではこの状態を「ばらつきが大きすぎる」過分散と呼んでいます。
sageへの入力:
# x=4での生存確率qを求める
q = _logistic(-2.15+0.51*4) print q }}
0.472527695655406 }}
sageへの入力:
# 図7.3(B)をプロット
lst_plt = list_plot(tbl4_lst) prd_plt = list_plot([_p(q, y, 8)*tbl4.sum() for y in (0..8)], plotjoined=True, rgbcolor="red") (lst_plt + prd_plt).show(figsize=5) }}
sageへの入力:
# x = 4 の平均と分散を求める(平均に対して分散が大きいことが分かる)
d4 = d[d.x==4] print d4.y.mean(), d4.y.var() }}
4.05 8.36578947368 }}
sageへの入力:
# 二項分布から期待されている分散は2なので、8.4はかなり大きい
N = 8 vr = N*q*(1-q) print vr }}
1.99396217995198 }}
個体差を考慮した分析手法として一般化線形混合分布(GLMM)が出てきます。
ここでは、個体\(i\)の個体差を\(r_i\)とし、以下の様なリンク関数を定義します。 $$ logit(q_i) = \beta_1 + \beta_2 x_i + r_i $$
ここで\(r_i\)をどのような分布にするかですが、久保本では平均0,分散\(s\)の正規分布と しています。(統計モデルとして便利だからと理由です)
glmmMLを使った分析では、個体毎に異なる独立パラメータをclusterオプションで指定します。 ここでは、idが個体を識別する指標として使われています。
sageへの入力:
r('library(glmmML)') }}
}}
sageへの入力:
# python版ではない、R一般線形混合モデルGLMMを使う
# 個体差は、平均0の分散sの正規分布に従うと仮定して先の問題を扱う
r('fit <- glmmML(cbind(y, N - y) ~ x, data = d, family = binomial, cluster = id)') }}
Call: glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = d, cluster = id)
(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
Residual deviance: 269.4 on 97 degrees of freedom AIC: 275.4 }}
GLMとの違いは、分析結果から予測される二項分布を見ると明らかです。 GLMのようなfitted.valuesがGLMMにはないので、分析結果の二項分布をmyfuncで定義して、 stat_function使ってプロットしてみました。
きれいに二項分布が表示されました。(これが図7.10の(A))
sageへの入力:
# 結果からロジスティック関数をmyfunc <- 1/(1 + exp(-(-4.190 + 1.005*z))*Nで定義
r('myfunc <- function(z) { 8*(1/(1 + exp(-(-4.190 + 1.005*z))))}')
# 結果をプロット
graph = preGraph("fig7.10-R.png") r('p <- ggplot(d, aes(x=x, y=y) )+ geom_point(position=position_jitter(width=.2, height=0), shape=21, size=2.5) + stat_function(fun=myfunc)') r('plot(p)') postGraph(graph) }}
もう一つGLMMから求められた\(r_i\)の正規分布をプロットすると以下の様になります。
sageへの入力:
# σ=2.4の正規分布
T = RealDistribution('gaussian', 2.4, seed=100) plot(lambda x : T.distribution_function(x), [x, -10, 10], figsize=5) }}
個体毎の尤度\(L_i\)は、以下の様に\(r_i\)を積分して求めることができます。 $$ L_i = \int_{-\infty}^{\infty} p(y_i | \beta_1, \beta_2, r_i) p(r_i | s) dr_i $$
GLMMの結果からqを求める関数_qを定義し、上記の積分をSageの数値積分を使って計算する関数_Liを以下の様に 作りました。rの積分範囲は、sの分布から-10から10としました。
チェックのために、久保本図7.8にあるx=4, r={-2.20, -0.60, 1.00, 2.60}のqの値と、 y={0, 4, 7}の尤度\(L_i\)を出力してみました。
最後に、Liに20(個体数)を掛けた分布をx=4のヒストグラムに重ね合わせてみました。 きれいに、図7.10(B)の通り、生存種子数を表現できているのが分かります。
このような積分を含む処理でも数式処理システムSageを使えば簡単に分析できるのが、 ご理解頂けたと思います。
sageへの入力:
# 回帰分析の結果からqを求める
def _q(x, r):
# 指定されたyの尤度を求める(上記分布から積分範囲は-10から10とした)
def _Li(y):
}}
sageへの入力:
_q(4, -2.2), _q(4, -0.6), _q(4, 1.0), _q(4, 2.6) }}
(0.0854891394348065, 0.316479106263684, 0.696354929823834, 0.919086532784535) }}
sageへの入力:
_Li(0), _Li(4), _Li(7) }}
(0.18814448426796998, 0.07913801041038665, 0.10840844640974827) }}
sageへの入力:
prd2_plt = list_plot([_Li(x)*20 for x in range(9)], plotjoined=True, rgbcolor="red") (lst_plt + prd2_plt).show(figsize=5) }}
皆様のご意見、ご希望をお待ちしております。