カウントデータは、平均が大きくなると分散も大きくなることが一般的に知られている。そのため、カウントデータをモデル化するときは、単回帰のように分散が常に一定の確率分布ではなく、平均が大きくなると分散も大きくなるような確率分布でモデル化する必要がある。ポアソン分布および負の二項分布などが、後者の分布にあたる。ポアソン分布は、平均と分散が同じ値になるのが特徴である。これに対して、負の二項分布は、分散が平均よりも大きな値をとることができるのが特徴である。そのため、ポアソン回帰で過分散問題が生じたとき、負の二項回帰でモデル化を試みることがある。このページでは、Stan を使用して、負の二項回帰でカウントデータをモデル化する例を示す。
サンプルデータとして、ガラパゴス島に生息している動物の種数データを使用する。このデータには、島の面積およびその島で生息している動物の種の数が記録されている。ここでは、島の面積とその島で生息している動物の種数を、負の二項回帰でモデル化する例を示す。
data(gala, package = 'faraway')
head(gala)
## Species Endemics Area Elevation Nearest Scruz Adjacent
## Baltra 58 23 25.09 346 0.6 0.6 1.84
## Bartolome 31 21 1.24 109 0.6 26.3 572.33
## Caldwell 3 3 0.21 114 2.8 58.7 0.78
## Champion 25 9 0.10 46 1.9 47.4 0.18
## Coamano 2 1 0.05 77 1.9 1.9 903.82
## Daphne.Major 18 11 0.34 119 8.0 8.0 1.84
plot(log10(gala$Area), gala$Species, xlab = 'log10(Area)', ylab = 'Species')
次に、Stan コードを使用して、負の二項回帰を記述していく。島の面積を x と種数を y とする。Stan の data ブロックで x および y を記述する。負の二項分布では、パラメーターとして μ（平均；mean）と φ（ばらつき；dispersion）の 2 つのパラメーターを受け取る。
平均 μ については、ポアソン回帰と同様に、線型予測子の値 βX = β0 + β1x を指数関数に代入することで計算する。そのため、Stan コードでは負の二項分布のパラメーター μ をそのままパラメーターとして定義するのではなく、transformed parameters ブロックで、μ = β0 + β1x のような変換を行なう。一方で、ばらつきは、そのままサンプリングによって求めるようにする。
入力データの定義、パラメーターの定義、パラメーター変換の定義のあとに、model ブロックでサンプリングによりパラメーター推定を行うモデル式を書く。また、パラメーターを推定した後に、島の面積が 1〜5 の範囲にあるときの種数の予測区間を計算したいので、generated quantities ブロックで再サンプリングを行うようにする。
data {
int N;
real x[N];
int y[N];
int new_N;
real new_x[new_N];
}
parameters {
// parameters for calculate mean of NB distribution
real beta_0;
real beta_1;
// dispersion of NB distribution
real phi;
}
transformed parameters {
real mu[N];
for (n in 1:N) {
mu[n] = exp(beta_0 + beta_1 * x[n]);
}
}
model {
for (n in 1:N) {
y[n] ~ neg_binomial_2(mu[n], phi);
}
}
generated quantities {
real yhat[new_N];
real muhat[new_N];
for (n in 1:new_N) {
muhat[n] = exp(beta_0 + beta_1 * new_x[n]);
yhat[n] = neg_binomial_2_rng(muhat[n], phi);
}
} 上の Stan コードで記述されたモデルは指数関数と負の二項分布の確率関数の両方を使用しているため、計算が不安定である。そのため、stan 関数を実行すると失敗する可能性がある。実行に失敗した場合は、何回がやり直せばそのうち 1 回は最後まで実行される。この不安定さを解消するために poisson_log 関数を使うことが推奨されている。
R を使用してデータを代入し負の二項回帰モデルのパラメーター推定を行う。また、島の面積の値はばらつきが大きいので、ここでは常用対数化を行って解析に用いる。ただし、そのまま対数化を行なうと、島の面積が負の値となるので、島の面積の単位を km2 から m2 に変換してから対数化を行なう。
library(rstan)
x <- log10(gala$Area * 1000000)
y <- gala$Species
new.x <- seq(1, 10, 0.1)
d <- list(x = x, y = y, N = length(x), new_x = new.x, new_N = length(new.x))
fit <- stan(file = 'nb.stan', data = d)
fit
## Inference for Stan model: nb.
## 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 -1.61 0.01 0.58 -2.69 -2.00 -1.62 -1.23 -0.41 1335 1
## beta_1 0.81 0.00 0.08 0.64 0.75 0.81 0.86 0.96 1733 1
## phi 2.62 0.02 0.74 1.45 2.10 2.52 3.05 4.25 1822 1
## mu[1] 79.59 0.18 10.72 60.90 71.99 78.95 86.06 103.31 3932 1
## mu[2] 27.73 0.07 3.85 21.24 25.09 27.39 30.09 35.81 3583 1
## ...
## muhat[89] 569.95 3.89 163.66 319.64 457.32 546.24 659.27 948.70 4320 1
## muhat[90] 619.19 4.36 182.79 341.50 493.53 592.81 719.27 1043.75 5402 1
## muhat[91] 672.74 4.89 204.08 364.87 532.03 643.71 782.40 1146.85 5103 1
## lp__ 10506.22 0.04 1.26 10502.88 10505.69 10506.53 10507.12 10507.62 1248 1
ms <- rstan::extract(fit, pars = 'yhat')
dim(ms$yhat)
## [1] 4000 91
df.pred <- data.frame(x = new.x,
lower = apply(ms$yhat, 2, quantile, prob = 0.025),
median = apply(ms$yhat, 2, quantile, prob = 0.500),
upper = apply(ms$yhat, 2, quantile, prob = 0.925))
g <- ggplot(df.pred, aes(x = x))
g <- g + geom_ribbon(aes(ymin = lower, ymax = upper), fill = '#000000', alpha = 0.4)
g <- g + geom_line(aes(y = median))
g <- g + geom_point(data = data.frame(x = x, y = y), aes(x = x, y = y))
g <- g + xlab('log10(Area)') + ylab('Species')
print(g)
95% 予測区間を確認すると、ポアソン回帰の結果とは異なり、大部分の観測点が予測区間の内側にあることが確認できる。このデータセットに関しては、ポアソン回帰よりも、負の二項回帰でモデル化した方が適しているといえる。