# Hierachical categoríes on a finite mixture model

**URL:** <https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907>\
**Category:** Modeling\
**Tags:** rstan, fitting-issues, specification\
**Created:** [September 8, 2020, 2:29am UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907 "2020-09-08T02:29:36Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![PavelGarcia](https://avatars.discourse-cdn.com/v4/letter/p/858c86/32.png) [@PavelGarcia](https://discourse.mc-stan.org/u/PavelGarcia)\
**Post date:** [September 8, 2020, 2:29am UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/1 "2020-09-08T02:29:36Z")

</div>

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?

```stan
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()

```

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [September 9, 2020, 5:11pm UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/2 "2020-09-09T17:11:05Z")

</div>

> [@PavelGarcia](#):
>
> `vector [n_groups] Theta; // Proportion of food sources on diet.`

What error are you getting? I assume first this should be:

```no-highlight
simplex[n_groups] theta; // Proportion of food sources on diet.

```

since `theta` needs to sum to 1.0 to work for a mixture.

---

<div class="post-metadata">

**Author:** ![PavelGarcia](https://avatars.discourse-cdn.com/v4/letter/p/858c86/32.png) [@PavelGarcia](https://discourse.mc-stan.org/u/PavelGarcia)\
**Post date:** [September 9, 2020, 6:23pm UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/3 "2020-09-09T18:23:44Z")

</div>

Thank you bbbales2.

I have implemented your suggestion before:

```stan
simplex[n_groups] Theta;

```

With it, I just can get the proportion for one specie, but I am looking to estimate Theta for 10 species at the same time.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [September 9, 2020, 8:14pm UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/4 "2020-09-09T20:14:00Z")

</div>

If you have one mixture model now, you could do:

```
simplex[n_groups] Theta[S]; // S being the number of species

```

And then you could fit `S` different mixture models, but my feeling is I’m missing the problem.

Can you describe in more detail what you are trying to do? Is it `S` different mixture models over `n_groups` elements? Or is it something else?

---

<div class="post-metadata">

**Author:** ![PavelGarcia](https://avatars.discourse-cdn.com/v4/letter/p/858c86/32.png) [@PavelGarcia](https://discourse.mc-stan.org/u/PavelGarcia)\
**Post date:** [September 9, 2020, 8:40pm UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/5 "2020-09-09T20:40:53Z")

</div>

I haven try that.

My goal is estimate how much leaf litter, biofilm and pure algae are eating more or less 10 species of aquatic insects on three types of streams. To estimate the what they are eating I am using stable isotope data. On the model I am trying to implement the following equations:

\delta y\_{i, t, s} \sim Normal (\sum\_{k,s} ^{t} P\_{k, t, s} \times (\delta X\_{k, s} + \Delta X\_{t} \* L\_{t}), \sigma \_{k,t} )

where \delta y\_{i, t, s} is the signature for \delta ^{13}C or \delta ^{15}N of each insect species (t), at each stream type (s); \delta X\_{k, s} is the stable isotope signature for the food sources (k) at each stream type (s); P\_{k, t} are the proportions of each type of food source (k) for each species (t) at each stream type (s). \Delta X\_t is a trophic enrichemnt factor for each stable isotope, and L\_{t} is the trophic level for each species (t).

So, as I understand I should combine a hierarchical model and a finite mixture model to do it. I have to start with a model where I have just a stream and just a species. Now, I am trying to integrate on the code the estimations of P for all the species that I am working on. Then I will try to incorporate the last level of stream types.

I have code n\_groups as the different elements on diet.  
I have ran your suggestion of

````stan
simplex [n_groups] Theta [Taxa_no] // where n_groups are the elements of diet to get proportions and Taxa_no the different species for I am estimating Theta. However, I got  
an error message on the model:

```stan
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
No matches for: 

  vector[] ~ dirichlet(vector)

Available argument signatures for dirichlet:

  vector ~ dirichlet(vector)

Real return type required for probability function.
 error in 'model2ba03e2f34f8_mod_exp5' at line 44, column 47
  -------------------------------------------------
    42: }
    43:   
    44: Theta ~ dirichlet(rep_vector(1.0, n_groups)); // Uniform prior
                                                      ^
    45:   
  -------------------------------------------------

````

---

<div class="post-metadata">

**Author:** ![PavelGarcia](https://avatars.discourse-cdn.com/v4/letter/p/858c86/32.png) [@PavelGarcia](https://discourse.mc-stan.org/u/PavelGarcia)\
**Post date:** [September 9, 2020, 10:35pm UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/6 "2020-09-09T22:35:17Z")

</div>

After following the advice of bbbales2,

The model is the following:

```stan

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
    
    simplex [n_groups] Theta [Taxa_no]; // 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[Taxa_no] ~ 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 (t in 1:Taxa_no){
                  for(k in 1:n_groups, Taxa_no) {
            contributions[t,k] = log(Theta[t,k]) + normal_lpdf(y_C[i] | (mu_C[k]+(Delta_C[t]*L[t], sigma_C[k]);
              }
          }
    target += log_sum_exp(contributions);
    }
  
  for(i in 1:N) {
        for (t in 1:Taxa_no){
              for(k in 1:n_groups) {
          contributions[t,k] = log(Theta[k]) + normal_lpdf(y_N[i] | (mu_N[k]+(Delta_N[t]*L[t])), sigma_N[k]);
                 }
           }
  target += log_sum_exp(contributions);
  }

  }

  
  "
,fill=TRUE)
sink()

```

However, it is not running. I am getting a syntax error of too many indexes on the likelihood part.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [September 10, 2020, 5:28pm UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/7 "2020-09-10T17:28:15Z")

</div>

> [@PavelGarcia](#):
>
> δyi,t,s∼Normal(∑tk,sPk,t,s×(δXk,s+ΔXt∗Lt),σk,t)

I don’t think the twiddle notation does very well with mixture models. Also I think the indices need a bit of work (\sum\_{k, s}^t seems confusing – for measurement y\_{i, t, s} sum to the t th species?).

Let’s also assume capital letters are the total number of things. So there are K food sources and T species and S stream types.

So let’s start with something like this:

p(y\_{i, t, s} | X, L, \sigma) = \sum\_k^K w\_{k, t, s} p\_k(y\_{i, t, s} | X\_k, L, \sigma\_k)\\

In this case I am assuming each component p\_k(y\_{i, t, s} | X\_k, L, \sigma\_k) is a normal likelihood with mean X\_{k, s} + \Delta X\_t L\_t and standard deviation \sigma\_{k, t}. The mixture weights w\_{k, t, s} are modeled hierarchically, and represent TS different simplexes of length K.

Anyway, start with that notation and push it in the direction you want it to go.

---

<div class="post-metadata">

**Author:** ![PavelGarcia](https://avatars.discourse-cdn.com/v4/letter/p/858c86/32.png) [@PavelGarcia](https://discourse.mc-stan.org/u/PavelGarcia)\
**Post date:** [September 10, 2020, 7:22pm UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/8 "2020-09-10T19:22:42Z")

</div>

Thank you bbbales2,

I will follow your instructions. I am wondering, if w\_{kt,s} is a matrix or a cov\_matrix?

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [September 10, 2020, 7:32pm UTC](https://discourse.mc-stan.org/t/hierachical-categories-on-a-finite-mixture-model/17907/9 "2020-09-10T19:32:09Z")

</div>

> [@PavelGarcia](#):
>
> is a matrix or a cov\_matrix

I’d represent it in Stan like:

```no-highlight
simplex[K] w[T, S];

```

That is a 2d array (dimension T by S) of K-dimensional simplexes.

Anyway I do think it’s valuable to try to write down your model in equations. Once you have everything clear there then it’s easier to translate to code. If you get stuck ask away.