Please share your Stan program and accompanying data if possible.

Hello all,

I have been trying to implement a hierarchical Dirichlet-Multinomial model for analyzing spatial variation in fish diet composition. My parameter of interest is the proportion of a given prey at a given site (theta in the model string). The softmax component in the hierarchical structure accounts for the influence of fish size on the probability of observing a prey. My model runs OK if I estimate the general diet composition across all locations. However, when I include the categorical covariate (site) in the model and do the necessary adaptations, Iâ€™m getting the following:

Rejecting initial value:

Chain 1: Gradient evaluated at the initial value is not finite.

Chain 1: Stan canâ€™t start sampling from this initial value.

Chain 1:

Chain 1: Initialization between (-2, 2) failed after 100 attempts.

Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

I have included the model specification, as well as some code to generate dummy data:

```
data {
int N_obs; // number of observations
int N_prey; // number of categories
int site[N_obs] ; // grouping factor
int N_sites; // number of sites
int y[N_obs]; // outcomes
int TL[N_obs];
}
parameters {
simplex[N_prey] theta[N_sites];
vector[N_prey] a;
vector[N_prey] b;
}
transformed parameters {
vector[N_prey] alpha;
vector[N_prey] s;
for (i in 1:N_obs) {
for (k in 1:N_prey) {
s[k] = a[k] + b[k]*TL[i];
alpha = softmax(s);
}
}
}
model {
a ~ normal(0,0.001);
b ~ normal(0,0.001);
for (k in 1:N_prey) {
theta[k] ~ dirichlet(alpha);
}
for (i in 1:N_obs) {
y[i] ~ categorical(theta[site[i]]);
}
}
```

# Simulating data and running the model

```
require(rstan)
require(rethinking)
N <-500
TL <-round(runif(N),1)*100
site <- c(rep(1,100), rep(2,100), rep(3,100), rep(4,100), rep(5,100))
b <-c(-0.5,0.4,1)
y <-rep(NA,N)
for (i in 1:N){
score <-0.5*(1:3)+b*TL[i]
p <-softmax(score[1],score[2],score[3])
y[i] <-sample(1:3,size=1,prob=p)
}
stan_data <- list(N_obs = N, N_prey = length(unique(y)), y = y,
TL = TL, site = site, N_sites = length(unique(site)))
stan.fit <- rstan::stan(file = "DMM.stan", data = stan_data,
warmup = 100,iter = 500, chains = 2)
```

Iâ€™d very much appreciate some help on that! Thanks in advance.

Edited by @jsocolar for R syntax highlighting.