I’m trying to estimate the effect of a predictor on count data. I also want to account for time-dependent trends in the counts by including 3 autoregressive terms. I input these t-k data instead of creating them in the stan model.

These data are grouped, so I want to incorporate varying parameter estimates and intercepts for each group.

When I run this model, I get many divergent transitions unless I run it for many iterations (~20,000). I think the likelihood surface is too complex and causing the “funnel” problem. A standard recommendation is to then implement non-centered parameterization. All the documentation/examples I can find for non-centered parameterization are for normally distributed data however, and I can’t quite figure out how to implement it for a negative binomial model.

Does this reparameterization even make sense? Or maybe there’s a better way to get this to fit more efficiently?

```
data {
int<lower=0> N; //number of observations
int count[N]; //Outcome variable
int Ngroups; //number of groups
int groupID[N]; //group ID
real predictor[N]; //Predictor variable
real ARcount1[N]; // count t-1
real ARcount2[N]; //count t-2
real ARcount3[N]; //count t-3
}
parameters {
real Intercept; //global intercept
real AR_term1[Ngroups]; //parameter for 1st autoregressive term
real AR_term2[Ngroups]; //parameter for 2ndst autoregressive term
real AR_term3[Ngroups]; //parameter for 3rdst autoregressive term
real<lower=0> sigma_gr; // sigma for estimating the group intercepts
real slope[Ngroups]; //the regression parameters on the predictor
real ran_intercept[Ngroups]; //Intercept for the groups
real<lower=0> phi; //the overdispersion parameters
}
model {
phi ~ cauchy(0, 2.5);
sigma_gr ~normal(0,1);
Intercept ~ normal(0,1); //priors following Gelman 2008
slope1 ~ cauchy(0,2.5);
AR_term1 ~ cauchy(0,2.5);
AR_term2 ~ cauchy(0,2.5);
AR_term3 ~ cauchy(0,2.5);
for (n in 1:N){
count[n] ~ neg_binomial_2(exp(Intercept + AR_term1[groupID[n]]*ARcount1[n]+
AR_term2[groupID[n]]*ARcount2[n]+AR_term3[groupID[n]]*ARcount3[n]+
slope1[groupID[n]]*predictor[n]+ran_intercept[groupID[n]]), phi);
}
for(j in 1:Ngroups){
ran_intercept[j]~normal(0,sigma_gr); //varying intercept by group
}
}
generated quantities { //generate within-sample predictions for validation
int<lower = 0> count_pred[N];
for (i in 1:N){
count_pred[i] = neg_binomial_2_rng(exp(Intercept + AR_term1[groupID[i]]*ARcount1[i]+
AR_term2[groupID[i]]*ARcount2[i]+AR_term3[groupID[i]]*ARcount3[i]+
slope1[groupID[i]]*predictor[i]+ran_intercept[groupID[i]]), phi);
}}
```