Hello Stan folks,

I’m trying to run the following hierarchical model with the attached stan model. However, I still encounter divergence problem after non-centered parameterization (I hope I did the non-centered parameterization right).

Would you please take a look on my code and possibly help me on this? Many thanks.

The model in my mind…

\beta_{zpSR} \sim b_0 + b_1\times Omni + b_2\times Omni^2

\mu_{phyDen} \sim X\times\beta + X_{zpSR} \times \beta_{zpSR}

phyDen \sim N(\mu_{phyDen}, \sigma_{phyDen})

The stan model I wrote…

```
data {
int<lower=1> N; // the number of observations
int<lower=1> K1; // number of column of the model matrix (5)
matrix[N, K1] X; // N (~1200 obs; ) by K1 (5 predictors) model matrix
vector[N] phyDen; // response variable, Y
vector[N] zpSR; //
int<lower=1> K2; // number of column of the model matrix determining the effects of zpSR (3)
matrix[N, K2] X2; // the model matrix determining the effects of zpSR
}
parameters {
vector[K1] gamma;
vector<lower=0>[K1] tau;
vector[K1] beta_raw[N];
vector[K2] gamma_beta_zpSR;
vector<lower=0>[K2] tau_beta_zpSR;
vector[K2] b_raw[N];
real<lower=0> sigma_phyDen; //standard deviation of the individual observations
}
transformed parameters{
vector[K1] beta[N];
vector[K2] b[N];
vector[N] mu_phyD;
vector[N] beta_zpSR;
for(i in 1:N){
beta[i] = gamma + tau .* beta_raw[i];
b[i] = gamma_beta_zpSR + tau_beta_zpSR .* b_raw[i];
}
for (i in 1:N){
beta_zpSR[i] = X2[i]*b[i];
mu_phyD[i] = X[i]*beta[i] + zpSR[i]*beta_zpSR[i];
}
}
model {
//priors
gamma ~ normal(0,5); //weakly informative priors on the regression coefficients
tau ~ cauchy(0,2.5); //weakly informative priors
gamma_beta_zpSR ~ normal(0,5); //weakly informative priors on the regression coefficients
tau_beta_zpSR ~ cauchy(0,2.5); //weakly informative priors
sigma_phyDen ~ gamma(2,0.1); //weakly informative priors
//likelihood
for (i in 1:N){
beta_raw[i] ~ normal(0, 1);
b_raw[i] ~ normal(0, 1);
}
phyDen ~ normal(mu_phyD, sigma_phyDen);
}
generated quantities {
vector[N] log_lik;
vector[N] y_pred;
for (n in 1:N){
log_lik[n] = normal_lpdf(phyDen[n] | mu_phyD[n], sigma_phyDen);
}
for (n in 1:N){
y_pred[n] = normal_rng(mu_phyD[n], sigma_phyDen);
}
}
```

The data is here if necessary…