My master's thesis involves modelling fMRI data. Each of $M$ participants has a total of $N$ voxels being measured. All these measurements represent an activation amplitude. Aside from this, I have data on participant characteristics (e.g. age, gender, etc.). My aim is to get posterior distributions across the assignment probabilities of the classes. In essence I want to cluster them using a mixture model. And I want to know how likely each cluster is per voxel. Either the voxel is negatively, positively, or 'null' activated, so there are three classes.
What I am trying to do is fit a three-component Gaussian mixture model. Let's denote this as follows using the allocation variable perspective, so I condition on some grouping $z_i$. For simplicity's sake I will just take the variance to be equal across components:
$$ y_i | z_i \sim \mathcal{N}(\mu_{z_i}, \sigma^2)\quad \text{with}\quad \pi(z_i = g) = \frac{1}{3} $$
As I understand each $\mu_{z_i}$ can again be modelled using a linear combintation of mixed effects. This is where I get stuck. As I undersand random-effects, they allow to account for group-level variation. For instance, each of my $M$ participants can be assigned its own random-intercept, allowing for prediction to be more accurate. However, this appears to be at odds with the formulation of the mixture model: there are only three components. How can I then get $M$ different values given the restriction of three components? What am I overlooking? My understanding of the mixture model is probably wrong.