【最尤法】ポアソン回帰の実装の仕方

マーケティングでポアソン回帰を使った事例があまり見つからなかったのですが、この本が分かりやすいかと思います。

マーケティングの統計モデル (統計解析スタンダード)

• 作者: 佐藤 忠彦
• 出版社/メーカー: 朝倉書店
• 発売日: 2015/08/25
• メディア: 単行本（ソフトカバー）

#######3章-ポアソン回帰モデルによる市場反応分析#######
#読み込みデータsec3_poisson_reg.csvのファイルレイアウト
#購買日付
#商品Aの販売個数
#商品Aの価格掛率
#商品Bの価格掛率
#商品Aの山積み陳列実施の有無（1実施0非実施）
#商品Bの山積み陳列実施の有無（1実施0非実施）
#来店客数
Dataset <- read.table("sec3_poisson_reg.csv",header=TRUE, sep=",", 
    na.strings="NA", dec=".", strip.white=TRUE)
##オフセット変数なし
GLM.2 <- glm(Sale_Unit_A ~ PriceIndex_A + PriceIndex_B + Display_A + 
    Display_B,family=poisson(log), data=Dataset)
summary(GLM.2)
##オフセット変数あり
GLM.3 <- glm(Sale_Unit_A ~ PriceIndex_A + PriceIndex_B + Display_A + 
    Display_B,family=poisson(log), offset=log(Visitors), data=Dataset)
summary(GLM.3)

結果は、下記となっています。

> summary(GLM.2)

Call:
glm(formula = Sale_Unit_A ~ PriceIndex_A + PriceIndex_B + Display_A + 
    Display_B, family = poisson(log), data = Dataset)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.1719   -1.3862   -0.1740    0.9973   19.8663  

Coefficients:
             Estimate Std. Error  z value Pr(>|z|)    
(Intercept)   6.35932    0.10355   61.412   <2e-16 ***
PriceIndex_A -6.96169    0.06148 -113.240   <2e-16 ***
PriceIndex_B  2.04395    0.08974   22.778   <2e-16 ***
Display_A     0.85543    0.01923   44.476   <2e-16 ***
Display_B    -0.07092    0.03073   -2.308    0.021 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 25654  on 1313  degrees of freedom
Residual deviance:  6188  on 1309  degrees of freedom
AIC: 10926

Number of Fisher Scoring iterations: 5

続いて、ポアソン回帰をもう少し深く理解するために、尤度関数を書いて、最尤法を使って解くことにします。

少しだけコードの説明。
ポアソン分布のパラメータはλですが、それは、exp(β0 + β1 * x1 + …)という線形の関係で書くことができます。
ポアソン分布の式から、尤度関数を求め、そのlogを取ったものが対数尤度関数になります。

log(y!)を計算するのですが、直接
log(gamma(y))
とすると、yが大きすぎて階乗の計算ができません。
大きいlog gammaの計算はlgammaで計算できます。

hh <- nrow(Dataset)  ## レコード数

## 尤度関数を用いた離散選択モデルを実施
# パラメータの初期化
beta_0 <- c(0, 0, 0, 0, 0)

# ポアソン回帰モデルの対数尤度関数
fr <- function(x)
{
  b0 <- x[1] #切片
  b1 <- x[2] #PriceIndex_A
  b2 <- x[3] #PriceIndex_B
  b3 <- x[4] #Display_A
  b4 <- x[5] #Display_B

  LL = 0

  for(i in 1:hh) {
    # 効用の計算
    lambda <- exp(b0 + b1 * Dataset[i, 3] + b2 * Dataset[i, 4] + 
        b3 * Dataset[i, 5] + b4 * Dataset[i, 6])

    # 対数尤度の計算
    LLL <- Dataset[i, 2] * log(lambda) -lambda -lgamma(Dataset[i, 2])
    LL <- LL + LLL
  }
  return(LL)
}

# 対数尤度関数の最大化
#  "fnscale = -1" で最大化
logit.out <- optim(beta_0, fr, method = "BFGS", hessian = TRUE, 
        control = list(fnscale = -1))

# 収束判定
logit.out$convergence

# 推定されたパラメータ
# パラメータの最尤推定値
logit.out$par

optimを使った最尤法の結果とglmを使ったパラメータの推定結果が一致することが確認できました。