Hello everyone,

I am trying to fit the model in the picture on simulated data. I used the same model for generating data except that the covariance matrix was an identity matrix, and the beta bars were standard normal random variables. I have tried both centered and non-centered parameterization; I encountered with following problems.

Firstly, there were some divergent samples in the centered model. Before the switch to the non-centered version, I tried the handle the problem by increasing the parameter of LKJ distribution and decreasing the step size. I worked and posterior beta distributions almost centered around the generated beta value. But Ä±t was not the same for beta bars. Their posterior distributions are far away from generated beta bars. I wonder how I can handle this problem.

Secondly, I tried non-centered parameterization. It worked, but run time was ten times higher than centered parameterization for the same number of iterations. Is that normal? I thought run time would decrease.

Also, chains did not mix, and almost all iterations exceeded the max_treedepth even for max_treedepth=13.

I shared the data simulation code and centered and non-centered models in this order. I would be very happy if anyone could help me.

### sim

```
data {
/* Dimensions */
int<lower=0> M; // number of varaibles
int<lower=0> N; // number of observations
/* Design Matrix */
matrix[N, M] X;
/* offset */
vector<lower=0>[N] offset;
/* Hyperparameters*/
vector<lower=0>[M] sigma; // diagonal values of the covariance matrix, equals to ones matrix.
}
generated quantities {
vector[M] betabar; // mean vector of multivariate normal distribution
vector[M] beta; // regression coefficients
real<lower=0> a; // prior for relative risk
vector<lower=0>[N] wi; //relative risk
vector[N] xi; // mean of poisson distribution
array[N] int<lower=0> y_rep; // response variable
for (i in 1 : M) {
betabar[i] = normal_rng(0, 1);
}
beta = multi_normal_rng(betabar, diag_matrix(sigma));
a = gamma_rng(1, 0.01);
for (j in 1 : N) {
wi[j] = gamma_rng(a, a);
}
xi = exp(X * beta + log(offset) + log(wi));
for (z in 1 : N) {
y_rep[z] = poisson_rng(xi[z]);
}
}
```

### NB2 centered

```
data {
/* Dimensions */
int<lower=0> M; // number of varaibles
int<lower=0> N; // number of observations
/* Design Matrix */
matrix[N, M] X;
/* offset */
vector<lower=0>[N] offset;
/* Outcome */
array[N] int<lower=0> y;
/* Hyperparameters*/
vector<lower=0>[M] sigma; //diagonal values of the covariance matrix, equals to 10s vector
}
parameters {
vector[M] beta; // coeffiencts in linear equation
real<lower=0> a; // prior for relative risk
vector<lower=0>[N] wi; // relative risk
vector[M] betabar; // mean vector of multivariate normal distribution
cholesky_factor_corr[M] Lcorr; // cholesky factor (L_u matrix for R)
}
transformed parameters {
vector[N] xi; // mean vector of poisson distribution
corr_matrix[M] R; // correlation matrix
cov_matrix[M] Sigma; // covariance matrix
xi = exp(X * betabar + log(offset) + log(wi));
R = multiply_lower_tri_self_transpose(Lcorr);
Sigma = quad_form_diag(R, sigma);
}
model {
target += normal_lpdf(betabar | 0, 1);
target += lkj_corr_cholesky_lpdf(Lcorr | 4);
target += multi_normal_lpdf(beta | betabar, Sigma);
target += gamma_lpdf(a | 1, 0.01);
target += gamma_lpdf(wi | a, a);
target += poisson_lpmf(y | xi);
}
```

### Non-centered

```
data {
/* Dimensions */
int<lower=0> M; // number of varaibles
int<lower=0> N; // number of observations
/* Design Matrix */
matrix[N, M] X;
/* offset */
vector<lower=0>[N] offset;
/* Outcome */
array[N] int<lower=0> y;
/* Hyperparameters*/
vector<lower=0>[M] sigma; //diagonal values of the covariance matrix, equals to 10s vector
}
parameters {
vector[M] delta; // std normal variable that multiplied with covariance matrix
vector[M] delta2; // std normal variable that multiplied with mean
real<lower=0> a;
vector<lower=0>[N] wi; // relative risk
cholesky_factor_corr[M] Lcorr; // cholesky factor (L_u matrix for R)
}
transformed parameters {
corr_matrix[M] R; // correlation matrix
vector[M] beta; // regression coefficients
vector[N] xi; // mean of poisson
R = multiply_lower_tri_self_transpose(Lcorr);
beta = delta2 * 10 + quad_form_diag(R, sigma) * delta;
xi = exp(X * beta + log(offset) + log(wi));
}
model {
target += normal_lpdf(delta2 | 0, 1);
target += lkj_corr_cholesky_lpdf(Lcorr | 2);
target += normal_lpdf(delta | 0, 1);
target += gamma_lpdf(a | 1, 0.01);
target += gamma_lpdf(wi | a, a);
target += poisson_lpmf(y | xi);
}
```