I’m in the process of migrating a JAGS code to STAN to compare performance. The lowest level of the model is Negative Binomial, but I have in JAGS used a transformation such that the distribution is defined by the mean and coefficient of variation of the equivalent Poisson-Gamma mixture. I’m getting confused about how the Negative Binomial is parameterized in STAN. Specifically, what should I replace the dnegbin expression with in the below code? (…which I appreciate may have other issues as I’m in the process of recoding it for STAN.)
data {
int<lower=1> NrOfLans;
int<lower=1> NrOfKretsar;
int<lower=1> Reports;
int<lower=0> K[Reports];
int<lower=0> A[Reports];
int<lower=1> KretsNrSeqID[Reports];
int<lower=1> LansSeqID[NrOfLans];
int<lower=1> LanSeqID[Reports];
int<lower=1> KretsarSeqID[NrOfKretsar];
int<lower=1> KretsListLanSeqID[NrOfKretsar];
}
parameters {
vector[NrOfLans] logL;
vector[NrOfKretsar] logK;
real logS;
real<lower=0> lansprec;
real<lower=0> kretsprecalpha;
real<lower=0> kretsprecmean;
real logcmean;
vector[NrOfLans] logckrets;
}
transformed parameters {
real<lower=0> c[NrOfLans];
real<lower=0> m[NrOfKretsar];
for (il in 1:NrOfLans){
c[il] = exp( logcmean + logckrets[il] );
}
for (ik in 1:NrOfKretsar){
m[ik] = exp( logL[KretsListLanSeqID[ik]] + logK[ik] );
}
}
model{
for (r in 1:Reports){
# Number of killed animals in report r is negative binomial distributed,
# but defined from the mean (m[KretsNr[r]]*A[r]) and coefficient of variation (c[Lan[r]]) of the gamma distribution of the equivalent gamma-poisson mixture
K[r] ~ dnegbin( 1 / ( 1 + m[KretsNrSeqID[r]] * A[r] * c[LanSeqID[r]]^2 ),c[LanSeqID[r]]^-2);
}
for(ik in 1:length(KretsarSeqID)){
logK[ik] ~ dnorm( 0 , kretsprec[KretsListLanSeqID[ik]] ); # logK[Kretsar[ik]] is the circuit random effect of circuit Kretsar[ik]
}
for (il in 1:length(LansSeqID)){
logL[il] ~ dnorm( logS , lansprec ); # logL[Lans[il]] is the county level random effect of county L
kretsprec[il] ~ dgamma( kretsprecalpha , kretsprecalpha / kretsprecmean); # kretsprec[Lans[il]] is the precission parameter for the hierarchical distribution of circuit effects with county Lans[il]
logckrets[il] ~ dnorm(0,logcprec); # county level effect for c
}
logS ~ dnorm(0,.001);
logcmean ~ dnorm(0,.01);
logcprec ~ dgamma(.01,.01);
kretsprecalpha ~ dgamma(.01,.01);
kretsprecmean ~ dgamma(.01,.01);
lansprec ~ dgamma(.01,.01);
}