ゼロ過剰負の二項モデルはカウントデータを扱う統計モデルの一種で、特に目的変数の側に0が多く含まれる場合に適している。
カウントデータを扱うモデルで最も有名なものはポアソン分布を利用したものであるが、ポアソン分布は性質上平均パラメタλが小さくなると分散も小さくなるので、データ全体の平均値が小さくなるとoverfittingが生じやすくなる。負の二項分布を用いたモデルはこの問題を部分的に回避できるが、データ中に0が多い場合には依然として問題が残る(南 et al. 2013）。
こうした際には、
1. 0でない値をとる確率についての回帰
2. 0でない値を取ったという条件のもとで、その数についての回帰
というような、階層的な回帰モデルを考えることで解決できる。
今回は柔軟な事前分布を設定するためStan言語を用いて尤度関数を記述したが、最尤推定で良いならRのパッケージ"pscl"で簡単に推定でき、パッケージ"MuMIn"によるモデル平均や事後モデル確率の計算にも対応する。ただし、正則化で扱うようなパラメータの多いモデルの場合にはモデル空間の大きさが現実的な時間で計算可能な範囲を超える。
[関連記事]
変数選択法まとめ1https://ameblo.jp/yusaku-ohkubo/entry-12267318650.html
変数選択法まとめ2https://ameblo.jp/yusaku-ohkubo/entry-12273211598.html
モデル平均法によるp-valueの危険性 https://ameblo.jp/yusaku-ohkubo/entry-12293763421.html
以下に、使用したソースコードを記載する。
======================================
data {
int<lower=1> N;
real x1[N];
real x2[N];
//real v3[N];
real x4[N];
real x5[N];
real x6[N];
real x7[N];
real x8[N];
real x9[N];
real x10[N];
real x11[N];
real x12[N];
real x13[N];
int<lower=0> y[N];
}
parameters{
vector[14] alpha;
vector[14] beta;
vector<lower=0>[14] lambda1;
vector<lower=0>[14] lambda2;
vector<lower=0>[14] gamma1;
vector<lower=0>[14] gamma2;
real<lower=0.00001> delta;
real<lower=0> tau1;
real<lower=0> tau2;
}
transformed parameters{
vector[N] theta;
vector[N] lambda;
for(n in 1:N){
theta[n]=inv_logit(x1[n]*alpha[1]+x2[n]*alpha[2]+pow(x2[n],2)*alpha[3]+x4[n]*alpha[4]+x5[n]*alpha[5]+x6[n]*alpha[6]+x7[n]*alpha[7]+x8[n]*alpha[8]+x9[n]*alpha[9]+x10[n]*alpha[10]+x11[n]*alpha[11]+x12[n]*alpha[12]+x13[n]*alpha[13]+alpha[14]);
lambda[n]=exp(x1[n]*beta[1]+x2[n]*beta[2]+pow(x2[n],2)*beta[3]+x4[n]*beta[4]+x5[n]*beta[5]+x6[n]*beta[6]+x7[n]*beta[7]+x8[n]*beta[8]+x9[n]*beta[9]+x10[n]*beta[10]+x11[n]*beta[11]+x12[n]*beta[12]+x13[n]*beta[13]+beta[14]);
}
}
model{
/// 以下で事前分布を設定する。ここではHorseshoe+を利用する。
/// Horseshoeであれば、gamma1,gamma2を削除して、lambda1, lambda2をcauchy(0,1)に変更する。
/// Laplace priorの場合、beta[]とalpha[]をそれぞれdouble_exponential(0, lambda1),double_exponential(0, lambda2)
/// pow(lambda1,2)~gamma(0.001,0.001); ow(lambda2,2)~gamma(0.001,0.001); に変更する
for(d in 1:14){
beta[d] ~ normal(0, lambda1[d] * tau1);
}
for(d in 1:14){
alpha[d] ~ normal(0, lambda2[d] * tau2);
}
lambda1 ~ cauchy(0, gamma1);
lambda2 ~ cauchy(0, gamma2);
gamma1~cauchy(0,1);
gamma2~cauchy(0,1);
tau1 ~ cauchy(0, 1);
tau2 ~ cauchy(0, 1);
///priorの設定、以上
///以下で尤度関数を記述する
for (n in 1:N) {
if (y[n] == 0) {
# Bernoulli(0|theta) + Bernoulli(1|theta) * Poisson(0|lambda)
target +=log_sum_exp(
bernoulli_lpmf(0| theta[n]),
bernoulli_lpmf(1| theta[n]) + neg_binomial_2_lpmf(0| lambda[n],delta)
);
} else {
# Bernoulli(1|theta) * Poisson(y|lambda)
target +=
log(theta[n]) + neg_binomial_2_lpmf(y[n]| lambda[n],delta
);
} }
}
///WAICの算出に必要な対数尤度を以下で計算する
generated quantities {
vector[N] log_likelihood;
int index;
//real y_pred;
for(n in 1:N){
if (y[n] == 0) {
# Bernoulli(0|theta) + Bernoulli(1|theta) * Poisson(0|lambda)
log_likelihood[n] =log_sum_exp(
bernoulli_lpmf(0| theta[n]),
bernoulli_lpmf(1| theta[n]) + neg_binomial_2_lpmf(0| lambda[n],delta)
);
} else {
# Bernoulli(1|theta) * Poisson(y|lambda)
log_likelihood[n] =(bernoulli_lpmf(1| theta[n]) + neg_binomial_2_lpmf(y[n]| lambda[n],delta));
}}
//y_pred = normal_rng(index == 1 ? mu : 0, 1);
}
参考文献
南美穂子. (2013). ゼロの多いデータの解析: 負の 2 項回帰モデルによる傾向の過大推定.