sage/データ解析のための統計モデリング入門勉強ノート第7章

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

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

ここで紹介した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です。以下に久保本から引用します。

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

↑

個体差の大きなデータ †

７章の例題は、個体$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]

↑

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

プロットデータを見ると、とても１００個のデータがあるようには見えません。これは、種子数・葉数が整数で同じ点に複数のデータが重なっているからです。

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)

↑

GLMを使った分析 †

このデータを６章と同じようにGLMで分析してみます。予測値は直線となっており、６章のようなロジスティック曲線にはなっていません。

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)

久保本に従い、$x_i = 4$のデータを抽出し、そのヒストグラムとGLMから推定された二項分布を表示してみましょう。

以下の様にPandasの抽出機能を使ってx=4のデータを取り出し、頻度を求めます。これにロジスティック_logisticと二項分布_pの２つの関数を定義します。

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$と比べて４倍ほど大きな値となっています。久保本ではこの状態を「ばらつきが大きすぎる」過分散と呼んでいます。

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への入力:

# 二項分布から期待されている分散は２なので、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$をどのような分布にするかですが、久保本では平均０，分散$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を使う
# 個体差は、平均０の分散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)

↑

$r_i$の分布 †

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

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

sageへの入力:

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

↑

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)

↑

コメント †

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

sage/データ解析のための統計モデリング入門勉強ノート第7章

最新の5件

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

前準備 †

個体差の大きなデータ †

データの特徴と可視化 †

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

GLMを使った分析 †

一般化線形混合分布（GLMM）を使った分析 †

結果の可視化 †

\(r_i\)の分布 †

GLMMから求められた尤度 †

コメント †