Hi, I am working on a model to estimate the proportion (Theta) of food types on freshwater insects using stable isotope data. I am interested in estimate the diet proportions for at least 10 different species for three types of streams. I have started with a finite mixture model. However, I am getting troubles with how I should write the code to estimate Theta for each species. Would someone give me any advice on how to code the species factor on my model?

```
sink("mod_exp5.stan")
cat("
data {
int N; // number of datapoints (bugs)
vector[N] y_C; // d13C signature for bugs
vector[N] y_N; // d15N signature for bugs
int n_groups; // number of food sources
vector [n_groups] Mu_C; // d13C mean value of each source
vector [n_groups] Mu_C_sd; // d13C sd value of each source
vector [n_groups] Mu_N; // d15N mean value of each source
vector [n_groups] Mu_N_sd; // d15N sd value of each source
int <lower=0> Taxa [N]; // Species of Bugs
int <lower=0> Taxa_no; // number of species
}
parameters {
ordered [n_groups] mu_C; // Estimation of mean d13C for the mix source food
ordered [n_groups] mu_N; // Estimation of mean d15N for the mix source food
vector [n_groups] Theta; // Proportion of food sources on diet.
vector <lower=0> [Taxa_no] Delta_C; // Trophic factor of enrichment for Carbon
vector <lower=0> [Taxa_no] Delta_N; // Trophic factor of enrichment for Nitrogen
vector <lower=0> [Taxa_no] L; // trophic level
vector<lower = 0> [n_groups] sigma_C; // Variation on the model for Carbon
vector<lower = 0> [n_groups] sigma_N; // Variation on the model for Nitrogen
}
model {
vector[n_groups] contributions;
// priors
for (k in 1: n_groups){
mu_C[k] ~ normal( Mu_C[k], Mu_C_sd[k]); //values from field Data samples of sources
}
for (k in 1: n_groups){
mu_N[k] ~ normal( Mu_N[k], Mu_N_sd[k]); //values from field Data samples of sources
}
Theta ~ dirichlet(rep_vector(1.0, n_groups)); // Uniform prior
Delta_C ~ normal (0.39, 1.14); // from Post (2002)
Delta_N ~ normal (3.4, 0.99); // from Post(2002)
//L ~ uniform (1,10);
sigma_C ~ normal(0, 100);
sigma_N ~ normal(0, 100);
// likelihood
for(i in 1:N) {
for(k in 1:n_groups, Taxa_no) {
contributions[k] = log(Theta[k]) + normal_lpdf(y_C[i] | (mu_C[k]+(Delta_C[Taxa[i]]*L[Taxa[i]])), sigma_C[k]);
}
target += log_sum_exp(contributions);
}
for(i in 1:N) {
for(k in 1:n_groups, Taxa_no) {
contributions[k] = log(Theta[k]) + normal_lpdf(y_N[i] | (mu_N[k]+(Delta_N[Taxa[i]]*L[Taxa[i]])), sigma_N[k]);
}
target += log_sum_exp(contributions);
}
}
"
,fill=TRUE)
sink()
```