FrontPage

2014/04/27からのアクセス回数 1631

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/44/

また、私の公開しているSageのサーバ(http://www15191ue.sakura.ne.jp:8000/)にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

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

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

前回に続き、 データ解析のための統計モデリング入門 (以下、久保本と書きます) の第7章の例題をSageを使って再現してみます。

久保本もようやくテーマの個体差が大きい場合の例題に入りました。今回のグラフは図7.10です。以下に久保本から引用します。

Fig7.10.png

数式処理システム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')

個体差の大きなデータ

7章の例題は、個体\(x_i\)で観測された以下のデータを含みます。

  • 調査種子数\(N_i = 8\)
  • 生存種子数\(y_i\)
  • 葉数\(x_i\)
  • 個体識別番号id

データの特徴と可視化

いつもの通りデータを読み込み、その内容、統計情報をチェックします。

sageへの入力:

# 7章のデータを読み込む
d = pd.read_csv(DATA+'data.csv')
d.head()
   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への入力:

d.describe()
         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]

同じ点に重なるデータの可視化

プロットデータを見ると、とても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.

fig-7.2.png

残念ながら、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)

fig7.2-R.png

GLMを使った分析

このデータを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:
    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 ' ' 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への入力:

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)

fig7.3-R.png

久保本に従い、\(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):
    return 1/(1 + exp(-z))

sageへの入力:

# 二項分布を定義
def _p(q, y, N):
    return binomial(N, y)*q^y*(1-q)^(N-y)

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)

sage0.png

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)を使った分析

個体差を考慮した分析手法として一般化線形混合分布(GLMM)が出てきます。

ここでは、個体\(i\)の個体差を\(r_i\)とし、以下の様なリンク関数を定義します。 $$ logit(q_i) = \beta_1 + \beta_2 x_i + r_i $$

ここで\(r_i\)をどのような分布にするかですが、久保本では平均0,分散\(s\)の正規分布と しています。(統計モデルとして便利だからと理由です)

glmmMLを使った分析では、個体毎に異なる独立パラメータをclusterオプションで指定します。 ここでは、idが個体を識別する指標として使われています。

sageへの入力:

#r("install.packages('glmmML')")
r('library(glmmML)')
 [1] "glmmML"    "jsonlite"  "ggplot2"   "stats"     "graphics"  "grDevices" "utils"     "datasets"
 [9] "methods"   "base"    

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)


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

fig7.10-R.png

\(r_i\)の分布

もう一つGLMMから求められた\(r_i\)の正規分布をプロットすると以下の様になります。

  • 10, 10ではほぼ0になっています。

sageへの入力:

# σ=2.4の正規分布
T = RealDistribution('gaussian', 2.4, seed=100)
plot(lambda x : T.distribution_function(x), [x, -10, 10], figsize=5)

sage1.png

GLMMから求められた尤度

個体毎の尤度\(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):
    return _logistic(-4.190 + 1.005*x + r)
# 指定されたyの尤度を求める(上記分布から積分範囲は-10から10とした)
def _Li(y):
    return numerical_integral(lambda r: _p(_q(4, r), y, 8)*T.distribution_function(r), -10, 10)[0]

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)

fig7.10_B.png

コメント

選択肢 投票
おもしろかった 0  
そうでもない 0  
わかりずらい 0  

皆様のご意見、ご希望をお待ちしております。


(Input image string)


添付ファイル: filesage1.png 237件 [詳細] filesage0.png 236件 [詳細] fileFig7.10.png 253件 [詳細] filefig7.10-R.png 239件 [詳細] filefig7.10_B.png 238件 [詳細] filefig7.3-R.png 224件 [詳細] filefig7.2-R.png 217件 [詳細] filefig-7.2.png 229件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-04-27 (日) 15:27:24 (1098d)
SmartDoc