Hi everyone,

I would like to know if I can implement this model in a single Stan code. The aim is to build a mixture where the components come from a hierarchical model and the weights from another model. The density should look like this:

f_{Y_{jm}}(y)=\sum_{i=1}^4\pi_{jm}(k=i) * f_{Y_{jm[k=i]}}\left(y/k=i\right)

Let’s say Y_{jm[k=i]} is one individual score in subject k, class j and school m. I have n=1,...,N scores, k=1,...,4 subjects and J=1,...,J classes, and the classes are in m=1,...,M schools.

Data is aggregate so the outcomes at student level Y_{jm[k=i]} are unobserved. I have the average for each subject within a school, plus the total number of students and exams for each subject.

The aim is to model the distribution of individual scores in class j and school m, with consistent partial pooling in a subject component towards the overall mean from all the schools for that subject, plus different weights given the proportion of exams allocated to each of them (again with partial pooling towards the overall mean of the weights).

I thought to recover the individual scores, imputing \sigma_{jm[k=i]} and a distribution of Y_{jm[k=i]} (both based on external research) and then simulate data. Finally, I’d fit the mixture model on the simulated dataset (or datasets) without the need of running two models.

Any help will be appreciated :)