トップ > R

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

R による XML の CSV 化

R Groovy XML

R を使って XML の内容 （特定の要素のみ） を CSV ファイルへ出力してみました。

  • R 3.0.1

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

CSV 化の対象 XML は下記で、VALUE 要素の属性とテキストノード値を CSV 出力する事にします。

対象 XML ファイル data.xml
<?xml version="1.0" encoding="UTF-8"?>
<GET_STATS_DATA>
  <STATISTICAL_DATA>
    <DATA_INF>
      <VALUE tab="20" cat01="10" cat02="0010" cat03="020" time="2013000000">15</VALUE>
      <VALUE tab="20" cat01="20" cat02="0010" cat03="020" time="2013000000">27</VALUE>
      <VALUE tab="20" cat01="30" cat02="0010" cat03="020" time="2013000000">36</VALUE>
      <VALUE tab="20" cat01="40" cat02="0010" cat03="020" time="2013000000">66</VALUE>
      <VALUE tab="20" cat01="50" cat02="0010" cat03="020" time="2013000000">47</VALUE>
    </DATA_INF>
  </STATISTICAL_DATA>
</GET_STATS_DATA>

かなり簡略化してますが、次世代統計利用システム API で取得したデータがこのような構造になっていました。

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

まずは R に XML パッケージをインストールしておきます。

packages.install("XML")

XML のパース関数

XML を DOM へパースする関数は下記のような種類があります。

一見どれを使えばよいか分からなくなってしまいそうですが、xmlInternalTreeParse・xmlNativeTreeParse・xmlParse は xmlTreeParse で useInternalNodes を TRUE に指定した場合と同じです。 （つまり xmlTreeParse が処理の本体で、残りは useInternalNodes のデフォルト値を TRUE に変えただけ）

関数 useInternalNodes パラメータのデフォルト値
xmlTreeParse FALSE
xmlInternalTreeParse TRUE
xmlNativeTreeParse TRUE
xmlParse TRUE

ここで、useInternalNodes を TRUE に指定すると xpathApply や getNodeSet のような XPath 式を使う関数を使える C レベルの XML ノードが戻り値で返ってきます。

ちなみに、SAX を処理する関数は xmlEventParse です。

1. データフレーム化して CSV 出力

まずは、下記のような手順でデータフレームを作って CSV 出力してみました。

  • (1) xmlParse で DOM へパース
  • (2) getNodeSet で XPath 式を使って VALUE 要素のリストを取得
  • (3) sapply で属性やテキストノード値を個々にベクトル化 （xmlGetAttr で属性値を取得、xmlValue でテキストノード値を取得）
  • (4) data.frame で (3) をデータフレーム化
  • (5) write.csv で (4) の内容を CSV ファイルへ出力

なお、テキストノード値は strtoi で integer へ変換しています。

sample.R
library(XML)

doc <- xmlParse("data.xml")
items <- getNodeSet(doc, "//VALUE")

tab <- sapply(items, function(x) xmlGetAttr(x, "tab"))
cat01 <- sapply(items, function(x) xmlGetAttr(x, "cat01"))
cat02 <- sapply(items, function(x) xmlGetAttr(x, "cat02"))
cat03 <- sapply(items, function(x) xmlGetAttr(x, "cat03"))
time <- sapply(items, function(x) xmlGetAttr(x, "time"))
value <- sapply(items, function(x) strtoi(xmlValue(x)))

df <- data.frame(tab, cat01, cat02, cat03, time, value, stringsAsFactors = FALSE)

write.csv(df, file = "data.csv", row.names = FALSE)

なお、write.csv で row.names = FALSE としていますが、row.names = TRUE のままでは行番号の列が出力されてしまうのでご注意ください。

実行例
> R CMD BATCH sample.R

出力結果は下記の通りです。 integer に変換した value 列は " で囲まれずに出力されています。

出力結果 data.csv
"tab","cat01","cat02","cat03","time","value"
"20","10","0010","020","2013000000",15
"20","20","0010","020","2013000000",27
"20","30","0010","020","2013000000",36
"20","40","0010","020","2013000000",66
"20","50","0010","020","2013000000",47

2. 行列を CSV 出力

次は、データフレームを使わずに sapply で行列を作って CSV 出力してみました。

なお、下記のように sapply した際の結果は tab や cat01 等が行となってしまうので（行と列が逆）、t で転置行列化して CSV 出力しています。

ちなみに、このサンプルでは strtoi は無意味なので使っていません。

sample2.R
library(XML)

doc <- xmlParse("data.xml")
items <- getNodeSet(doc, "//VALUE")

d <- sapply(items, function(x) list(
    tab = xmlGetAttr(x, "tab"),
    cat01 = xmlGetAttr(x, "cat01"),
    cat02 = xmlGetAttr(x, "cat02"),
    cat03 = xmlGetAttr(x, "cat03"),
    time = xmlGetAttr(x, "time"),
    value = xmlValue(x)
))

write.csv(t(d), file = "data2.csv", row.names = FALSE)
実行例
> R CMD BATCH sample2.R

出力結果は下記の通りです。

出力結果 data2.csv
"tab","cat01","cat02","cat03","time","value"
20,10,0010,020,2013000000,15
20,20,0010,020,2013000000,27
20,30,0010,020,2013000000,36
20,40,0010,020,2013000000,66
20,50,0010,020,2013000000,47

3. SAX で CSV 出力

最後に SAX で CSV 出力してみました。 VALUE 要素のテキストノードを処理するために .state の値を制御フラグとして使っています。

また、write.csv は使わずに cat を使ってファイル出力しました。

sample3.R
library(XML)

f <- file("data3.csv", "w")

cat('"tab","cat01","cat02","cat03","time","value"\n', file = f, append = TRUE)

# VALUE 要素の属性を処理
procValueNode <- function(name, attrs, .state) {
    if (name == "VALUE") {
        cat(attrs[1], attrs[2], attrs[3], attrs[4], attrs[5], "", file = f, sep = ",", append = TRUE)
        # procValueText で処理させるために .state を TRUE に変更
        .state = TRUE
    }
    .state
}

# VALUE 要素のテキストノードを処理
procValueText <- function(content, .state) {
    if (.state) {
        cat(content, "\n", file = f, sep = "", append = TRUE)
    }
    # .state を FALSE とするために FALSE を返す
    FALSE
}

xmlEventParse("data.xml", handlers = list(
    startElement = procValueNode,
    text = procValueText
), state = FALSE)

close(f)
実行例
> R CMD BATCH sample3.R

出力結果は 2. と同じです。

出力結果 data3.csv
"tab","cat01","cat02","cat03","time","value"
20,10,0010,020,2013000000,15
20,20,0010,020,2013000000,27
20,30,0010,020,2013000000,36
20,40,0010,020,2013000000,66
20,50,0010,020,2013000000,47

処理時間の比較

10万行の CSV を出力するようなデータを処理した場合の処理時間は概ね下記のようになりました。比較のために Groovy で XmlSlurper や StAX を使った場合も試してみました。

タイプ スクリプト 処理時間
1. データフレーム化して csv 出力 sample.R 40秒
2. 行列を CSV 出力 sample2.R 40秒
3. SAX で CSV 出力 sample3.R 14秒
Groovy で XmlSlurper 利用 sample_xmlslurper.groovy 12秒
Groovy で StAX 利用 sample_stax.groovy 6秒

R の中ではやはり SAX を用いた方が高速でしたが、普通に Groovy で処理した方が速かったので、今回のようなケースを R で処理する必然性は低いかもしれません。

なお、使用したGroovy スクリプトは下記の通りです。

sample_xmlslurper.groovy
def doc = new XmlSlurper().parse(new File(args[0]))

println '"tab","cat01","cat02","cat03","time","value"'

doc.STATISTICAL_DATA.DATA_INF.VALUE.each {
    println "${it.@tab},${it.@cat01},${it.@cat02},${it.@cat03},${it.@time},${it.text()}"
}
sample_stax.groovy
import javax.xml.stream.*

def factory = XMLInputFactory.newInstance()
def xr = factory.createXMLStreamReader(new File(args[0]).newReader("UTF-8"))

def procValueNode = { stream ->
    if (stream.name.localPart == 'VALUE') {
        def items = (0..<stream.attributeCount).collect {
            stream.getAttributeValue(it)
        }
        items << stream.elementText

        println items.join(',')
    }
}

println '"tab","cat01","cat02","cat03","time","value"'

while(xr.hasNext()) {
    switch (xr.eventType) {
        case XMLStreamConstants.START_ELEMENT:
            procValueNode(xr)
            break
    }
    xr.next()
}

xr.close()