2015-02-23

以前、「R でポアソン回帰 - glm, MCMCpack」にて試した GLM によるポアソン回帰を Julia で実施してみました。

なお、Julia は開発中の v0.4.0 を使用しました。

Julia v0.4.0-dev

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20150223/

はじめに

今回は下記のパッケージを使用しますので、予めインストールしておきます。

パッケージ名	用途
DataFrames	csv ファイルの読み込み。データフレーム化
GLM	GLM の実施
Gadfly	グラフの描画
Cairo	PNGファイルへのグラフ描画

パッケージのインストール例

> julia
・・・

julia> Pkg.add("DataFrames")
・・・
julia> Pkg.add("GLM")
・・・
julia> Pkg.add("Gadfly")
・・・
julia> Pkg.add("Cairo")
・・・

ポアソン回帰を試すデータは以前に使ったものをそのまま使用します。

データ data3a.csv

y,x,f
6,8.31,C
6,9.44,C
6,9.5,C
12,9.07,C
10,10.16,C
・・・

データ内容は下記の通りです。

今回は、体サイズ x が種子数 y にどのように影響しているかをポアソン回帰で解析します。

項目	内容
y	種子数
x	植物の体サイズ
f	施肥処理

GLM によるポアソン回帰

y ~ x のモデルを使ってポアソン回帰を実施します。

以前、R で実施したものと同じポアソン回帰を行うには、glm 関数へ Poisson() と LogLink() （リンク関数）を指定します。

poissonGlm.jl

using DataFrames, GLM

d = readtable("data3a.csv")

res = glm(y ~ x, d, Poisson(), LogLink())

println(res)

実行結果は、以下のように R の結果とほぼ同じになりました。（前回の結果を参照）

今回は Julia の開発バージョンを使ったため deprecated syntax の WARNING がいくつか出力されました。

実行結果

> julia poissonGlm.jl
・・・
Coefficients:
              Estimate Std.Error z value Pr(>|z|)
(Intercept)    1.29172  0.363686 3.55174   0.0004
x            0.0756619 0.0356042 2.12509   0.0336

この結果より、平均種子数の予測値を求める関数は ${ \lambda = \exp(1.29172 + 0.0756619 x) }$ となります。

予測線の描画1

次に、ポアソン回帰の結果を使って種子数 y の予測線を描画し PNG ファイルへ出力してみます。

まず、種子数 y の予測値を算出するために、data3a.csv における x の最小値 7.19 から最大値 12.4 までを 0.1 刻みにした配列 xx を作成しました。

coef 関数を使えば glm の結果から Coefficients の Estimate 値を配列として取得できるので、これと xx を使って種子数 y の予測値を算出できます。

plot と layer を組み合わせる事で複数のグラフを重ねて描画できます。

poissonGlm_draw1.jl

using DataFrames, GLM, Gadfly

d = readtable("data3a.csv")

res = glm(y ~ x, d, Poisson(), LogLink())

# x の最小値 7.19 ～ 最大値 12.4 まで 0.1 刻みのデータを用意
xx = [minimum(d[:x]):0.1:maximum(d[:x])]
# 種子数 y の予測値を算出
yy = exp(coef(res)[1] + coef(res)[2] * xx)

# グラフの描画
p = plot(
    # data3a.csv の内容を描画
    layer(d, x = "x", y = "y", color = "f", Geom.point),
    # 予測線の描画
    layer(DataFrame(x = xx, y = yy), x = "x", y = "y", Geom.line)
)

# PNG 形式で出力
draw(PNG("poissonGlm_draw1.png", 500px, 400px), p)

f:id:fits:20150223203652p:plain

ここで、xx は以下のような内容です。

xx の内容

53-element Array{Float64,1}:
  7.19
  7.29
  ・・・
 12.29
 12.39

予測線の描画2

最後に、種子数 y の予測値を predict 関数で算出し、同様のグラフを描画します。

predict を使うには x の値だけでは無く (Intercept) の Estimate 値へ乗算する値も必要になるので、下記では DataFrame （下記の nd）へ項目 n （常に 1）として用意しました。

predict 内では ${ \exp(1.29172 n + 0.0756619 x) }$ のような計算が実施されます。

poissonGlm_draw2.jl

・・・
res = glm(y ~ x, d, Poisson(), LogLink())

# 7.19 ～ 12.4 まで 0.1 刻みのデータを用意
xx = [minimum(d[:x]):0.1:maximum(d[:x])]
# predict へ渡すデータを作成
nd = DataFrame(n = [1 for i = 1:length(xx)], x = xx)

# predict で種子数 y の予測値を算出
yy = predict(res, nd)

p = plot(
    # data3a.csv の内容を描画
    layer(d, x = "x", y = "y", color = "f", Geom.point),
    # 予測線の描画
    layer(DataFrame(x = xx, y = yy), x = "x", y = "y", Geom.line)
)

draw(PNG("poissonGlm_draw2.png", 500px, 400px), p)

f:id:fits:20150223203711p:plain

なんとなくな Developer のメモ

Julia でポアソン回帰 - glm