読者です 読者をやめる 読者になる 読者になる

R でポアソン回帰 - glm, MCMCpack

R GLM Bayes

書籍 「 データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学) 」の 3章 「一般化線形モデル(GLM)」 と 9章 「GLMのベイズモデル化と事後分布の推定」 で説明されていたポアソン回帰を下記のような 3通りで試してみました。

書籍では、R から WinBUGS を呼び出して MCMC サンプリングを行っていましたが、今回は R 上でベイズ統計解析を実施する MCMCpack パッケージを試してみました。

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

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

はじめに

まず、MCMCpack パッケージを使用するためにパッケージを R の実行環境へインストールしておきます。

MCMCpack のインストール
> install.packages("MCMCpack")

次に、ポアソン回帰を試すデータを書籍のサポート web サイトから取得してきます。

データ 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 と肥料による施肥処理 f が種子数 y にどのように影響しているかをポアソン回帰で解析します。

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

(1) GLM によるポアソン回帰 (glm 関数)

それでは glm 関数を使ったポアソン回帰を試してみます。 ここでは stepAIC() を使った AIC によるモデル選択を試してみました。

poissonGlm.R
d <- read.csv('data3a.csv')

# 説明変数を全て投入したモデル (y ~ x + f と同じ)
d.all <- glm(y ~ ., data = d, family = poisson)

library(MASS)
# AIC によるモデル選択
d.res <- stepAIC(d.all)

summary(d.res)

png("poissonGlm.png")
plot(d$x, d$y, col = c("red", "blue")[d$f])
xx <- seq(min(d$x), max(d$x), length = 1000)
# y ~ x のモデルを使った平均種子数の予測値の曲線を描画
lines(xx, exp(d.res$coefficients["(Intercept)"] + d.res$coefficients["x"] * xx), col="green")
# 以下でも可
#lines(xx, predict(d.res, newdata = data.frame(x = xx), type = "response"), col = "green")
dev.off()

なお、今回は y ~ x のモデルが選択される事 (施肥処理は効果が無い) を予め分かっているので、y ~ x のモデルを使って平均種子数の予測値を描画 (緑の線) しています。

ここで、 y ~ x モデルのリンク関数 (対数リンク関数) は { \log \lambda_i = \beta_1 + \beta_2 x_i } で、平均種子数 { \lambda_i }{ \lambda_i = \exp (\beta_1 + \beta_2 x_i) } となります。

また、{ \beta_1 }最尤推定値が d.res$coefficients["(Intercept)"]{ \beta_2 }最尤推定値が d.res$coefficients["x"] に該当します。

実行

実行すると下記のような結果になります。

> R CMD BATCH poissonGlm.R
実行結果 poissonGlm.Rout
・・・
> d.res <- stepAIC(d.all)
Start:  AIC=476.59
y ~ x + f

       Df Deviance    AIC
- f     1   84.993 474.77
<none>      84.808 476.59
- x     1   89.475 479.25

Step:  AIC=474.77
y ~ x

       Df Deviance    AIC
<none>      84.993 474.77
- x     1   89.507 477.29
> 
> summary(d.res)

Call:
glm(formula = y ~ x, family = poisson, data = d)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3679  -0.7348  -0.1775   0.6987   2.3760  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.29172    0.36369   3.552 0.000383 ***
x            0.07566    0.03560   2.125 0.033580 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 89.507  on 99  degrees of freedom
Residual deviance: 84.993  on 98  degrees of freedom
AIC: 474.77

Number of Fisher Scoring iterations: 4

・・・

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

f:id:fits:20131215205033p:plain

(2) MCMCpack を使ったベイズ統計によるポアソン回帰1 (MCMCpoisson 関数)

MCMCpack でポアソン回帰を行うには MCMCpoisson 関数を使うのが簡単だと思います。 基本的に glm 関数と同様のモデル式とデータを指定するだけです。

ここでは y ~ x のモデルだけをポアソン回帰してみました。

なお、glm 関数では線形予測子のパラメータ { \beta_1 }{ \beta_2 }最尤推定値を取得しますが、MCMCpoisson 関数では { \beta_1 }{ \beta_2 } のそれぞれの分布を取得する点が大きく異なります。

poissonMcmcPoisson.R
library(MCMCpack)

d <- read.csv('data3a.csv')

d.res <- MCMCpoisson(y ~ x, data = d)

summary(d.res)

png("poissonMcmcPoisson.png")
plot(d$x, d$y, col = c("red", "blue")[d$f])
xx <- seq(min(d$x), max(d$x), length = 10000)
lines(xx, exp(mean(d.res[,1]) + mean(d.res[,2]) * xx), col="green")
dev.off()

glm の時とは異なり、{ \beta_1 }{ \beta_2 } のそれぞれの算術平均値 (mean の結果) を使って、平均種子数の予測値を描画 (緑の線) しています。

実行

実行すると下記のような結果になります。

> R CMD BATCH poissonMcmcPoisson.R
実行結果 poissonMcmcPoisson.Rout
・・・
> summary(d.res)

Iterations = 1001:11000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean     SD Naive SE Time-series SE
(Intercept) 1.29315 0.3602 0.003602       0.010757
x           0.07547 0.0352 0.000352       0.001047

2. Quantiles for each variable:

                2.5%     25%    50%     75%  97.5%
(Intercept) 0.586606 1.05208 1.2956 1.54153 1.9917
x           0.006025 0.05125 0.0756 0.09951 0.1438

・・・

Mean の値が glm の結果とほぼ同じ値になっています。

f:id:fits:20131215205111p:plain

ちなみに、d.res には MCMC サンプリング結果として 1万個の { \beta_1 }{ \beta_2 } の値が格納されています。

{ \beta_1 } の分布

f:id:fits:20131215205124p:plain

{ \beta_2 } の分布

f:id:fits:20131215210052p:plain

(3) MCMCpack を使ったベイズ統計によるポアソン回帰2 (MCMCmetrop1R 関数)

最後に MCMCmetrop1R 関数を使ってみます。

MCMCmetrop1R 関数では自前で定義した尤度関数を使って MCMC サンプリングを実施できるので汎用的に使えます。

今回のモデルの対数尤度は { \log L ( \beta_1, \beta_2 ) = \sum_i \log \frac{\lambda_i ^ y_i \exp(- \lambda_i)}{y_i!} } で、{ \lambda_i = \exp (\beta_1 + \beta_2 x_i) } なので、これらを関数化 (下記の func) しています。

尤度関数が対数尤度か否かは logfun で指定するようになっており (TRUE の場合が対数尤度)、デフォルト値は TRUE となっています。 (今回は対数尤度なので logfun = TRUE です)

また、theta.init{ \beta_1 }{ \beta_2 } の初期値を c(0, 0) と指定しています。

poissonMcmcMetrop.R
library(MCMCpack)

d <- read.csv('data3a.csv')

# 尤度関数(対数尤度関数)
func <- function(beta, x, y) {
  lambda <- exp(beta[1] + beta[2] * x)
  sum(log(dpois(y, lambda)))
}

d.res <- MCMCmetrop1R(func, theta.init = c(0, 0), x = d$x, y = d$y, logfun = TRUE)
# 下記でも同じ
#d.res <- MCMCmetrop1R(func, theta.init = c(0, 0), x = d$x, y = d$y)

summary(d.res)

png("poissonMcmcMetrop.png")
plot(d$x, d$y, col = c("red", "blue")[d$f])
xx <- seq(min(d$x), max(d$x), length = 10000)
lines(xx, exp(mean(d.res[,1]) + mean(d.res[,2]) * xx), col="green")
dev.off()

実行

実行すると下記のような結果になります。

> R CMD BATCH poissonMcmcMetrop.R
実行結果 poissonMcmcMetrop.Rout
・・・
> summary(d.res)

Iterations = 501:20500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean      SD  Naive SE Time-series SE
[1,] 1.29545 0.35797 0.0025313      0.0077477
[2,] 0.07528 0.03498 0.0002473      0.0007561

2. Quantiles for each variable:

         2.5%     25%     50%     75%  97.5%
var1 0.584646 1.05368 1.29884 1.53850 1.9830
var2 0.007798 0.05159 0.07518 0.09907 0.1448

・・・

glm や MCMCpoisson と似たような結果となりました。

f:id:fits:20131215205205p:plain

なお、MCMCpoisson と MCMCmetrop1R ではバーンイン burnin とサンプリング数 mcmc のデフォルト値が異なっています。

関数 burnin のデフォルト値 mcmc のデフォルト値
MCMCpoisson 1000 10000
MCMCmetrop1R 500 20000