ゼロ過剰負の二項回帰のStanによる実装

今回は、ゼロ過剰負の二項回帰モデルをやりたいと思います。
なんか名前がゴツゴツしいですね。

とりあえず、今回もStanコードを書くのがメインなので、モデル部分はザックリいきます。

モデル
Stanによる実装

モデル

ゼロ過剰とは

ゼロ過剰ポワソン分布は有名ですね。
要は確率 $\omega$ でゼロしか出ない分布、確率 $1-\omega$ でポワソン分布が選ばれる混合分布を考える時ゼロ過剰ポワソン分布と呼びます。本来のポワソン分布よりも $\omega$ だけゼロが過剰に出ていますよね。だからゼロ過剰です。

では、ゼロ過剰負の二項分布とは何かというと、先ほどのポワソンを負の二項分布に変えればゼロ過剰負の二項分布です。

負の二項分布

今回は次のような定義で負の二項分布を定義しておきます。

$f_{NB} = \frac{\Gamma(\theta+Y_i)}{\Gamma(\theta)\Gamma(Y_i+1)}(\frac{\theta}{\theta+\mu_i})^{\theta} (\frac{\mu_i}{\theta+\mu_i})^{y_i}$

負の二項回帰

更に次のような構造を仮定することで回帰することにします。

$log(\mu_i) = \beta_0+\beta_1 X_{(1,i)} +\cdots + \beta _{Kb} X_{(Kb,i)}$
$log(\frac{\omega_i}{1-\omega_i}) =logit(\omega_i) =\gamma_0 + \gamma_1 x_{(1,i)} +\cdots + \gamma_{Kg} x_{(kg,i)}$

$\beta,\gamma$ がパラメータになります。

Stanによる実装

というわけでコードを書いてみました。

例の有名な教科書のゼロ過剰ポワソンの部分のコードを真似してます。
<br />

model="
functions{
	real ZINB_lpmf(int Y,real omega,real theta,real mu){
		if(Y==0){
			return log_sum_exp(bernoulli_lpmf(0|omega),bernoulli_lpmf(1|omega)*neg_binomial_2_lpmf(0|mu,theta));
		}else{
			return bernoulli_lpmf(1|omega)*neg_binomial_2_lpmf(Y|mu,theta);
		}
	}
}

data{
	int N;
	int Kb;
	int Kg;
	int<lower=0> Y[N];
	matrix[N,Kb] X;
	matrix[N,Kg] XX;
}

parameters{
	real<lower=0> theta;
	real beta0;
	vector[Kb] beta;
	real gamma0;
	vector[Kg] gamma;
}

transformed parameters{
	real<lower=0> mu[N];
	real<lower=0,upper=1> omega[N];
	real q[N];



	for(i in 1:N){
		q[i] = gamma0+XX[i]*gamma;
		omega[i] = inv_logit(q[i]);
		mu[i] =exp(beta0+X[i]*beta);

	}	
}

model{
	for(i in 1:N){
		Y[i] ~ ZINB(omega[i],theta,mu[i]);
	}
}"

収束しない場合は、弱事前分布ですかね。
回帰係数の部分に正規分布かt分布をあてるのが良いかと。