Stan によるポアソン回帰のパラメーター推定

ポアソン回帰

Stan によるポアソン回帰のパラメーター推定を行うために、ポアソン分布に従うデータを乱数生成する。ここで、小麦の個体 i の乾燥重量を xi とし、その個体の穂数を yi とし、両者は yi = Poison(λ = 2xi + 5) の関係を満たすように乱数生成する。

x <- rnorm(100, mean = 20, sd = 4)
y <- rpois(100, lambda = 2 * x + 5)
head(cbind(x, y))
##             x  y
## [1,] 10.50902 26
## [2,] 26.38858 51
## [3,] 26.58175 61
## [4,] 22.74905 53
## [5,] 17.33573 45
## [6,] 18.59018 40

GLM

R の glm 関数でポアソン回帰を行うときは、次のようにする。乱数データの生成条件とは異なるパラメーターが推測されたものの、図示すると、だいたいあっていることがわかる。

m <- glm(y ~ x, family = poisson(link = "log"))
m
## Call:  glm(formula = y ~ x, family = poisson(link = "log"))
## 
## Coefficients:
## (Intercept)            x
##     2.94791      0.04281
## 
## Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
## Null Deviance:	    259.3
## Residual Deviance: 84.77 	AIC: 650.4

plot(x, y)
lines(seq(min(x), max(x), 0.1), exp(coef(m)[[1]] + coef(m)[[2]] * seq(min(x), max(x), 0.1)))
ポアソン回帰の結果

Stan

Stan を利用してポアソン回帰のパラメーターを計算する。Stan のコードでは、まず個体の乾燥重量 x と穂数 y のデータを受け取るように data ブロックで定義する。次に、ポアソン回帰のパラメーターとして β0 および β1parameters ブロックで定義する。ポアソン回帰のモデル式では、直接に β0 および β1 を用いるわけではないので、β0 および β1 を λ に変換するように transformed parameters で指示する。最後に、model ブロックでサンプリングによりパラメーター推定を行うモデル式を書く。

data {
    int N;
    real x[N];
    int y[N];
}
parameters {
    real beta_0;
    real beta_1;
}
transformed parameters {
    real lambda[N];
    for (n in 1:N)
        lambda[n] = exp(beta_0 + beta_1 * x[n]);
}
model {
    y ~ poisson(lambda);
}

R のコードを以下のようにして、ポアソン回帰のパラメーター推定を行う。入力データの数分だけの λ も計算されるが、無視して構わない。Stan を利用したパラメーター推定では、β0 = 2.93、β1 = 0.04 と推定され、GLM と似た結果が得られた。

library(rstan)
d <- list(x = x, y = y, N = length(x))
fit <- stan(file = 'pois.stan', data = d)
fit
## Inference for Stan model: pois.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##                 mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
## beta_0          2.93    0.00 0.08     2.78     2.87     2.93     2.98     3.08   660    1
## beta_1          0.04    0.00 0.00     0.04     0.04     0.04     0.04     0.05   653    1
## lambda[1]      46.66    0.01 0.69    45.30    46.20    46.66    47.13    47.98  3919    1
## lambda[2]      61.00    0.06 1.78    57.55    59.78    60.96    62.18    64.44   921    1
## lambda[3]      51.78    0.02 0.95    49.95    51.14    51.78    52.41    53.61  1650    1