# Can anybody help with multi-level mixture models

Hello,

I am trying to fit a multi-level mixture model and running into a lot of problems. In short: I want to model the likelihood of participants making movements in various directions. I am trying to model this as a mixture of eight von Mises distributions (with fixed mu, and for now, fixed kappa).

When I try to extend this to a multi-level framework, I encounter a series of problems. First, I have realised that my standard method for specifying the correlation doesn’t work (as I have 8 directions x 6 conditions, so a 48x48 varcov matrix. For now, I have decided to remove this part of them model and instead just try and estimate the standard deviation of how each theta varies in each of my conditions.

I think my main source of issues is related to combining the fixed-effect theta with the random effects. The end result needs to be a simplex and I’m stumped on how to do this. Are there any standard tricks that I should know about?

My latest attempt involves defining theta as a simplex[8], and u[8] ~ normal(0, sd). Then I define u_theta = theta * exp(u). This results in theta>0 and it can lead to quite nice prior predictions, but the model doesn’t fit to data at all well. I suspect it may be due to sum(-theta) != 1 ? I could normalise all the u_theta, but I suspect this might lead to poor specification for my u parameters.

I’m quite close to giving up on the multi-level implementation and simply fitting the model to each person one at a time.

Code below… I hope it makes sense. There’s a good chance I’ve missed something daft.

Thanks in advance… I really appreciate how helpful people are on this forum.

``````data {
int<lower = 0> N;
int <lower = 1> K; // number of experimental conditions
array[N] int <lower = 1, upper = K> X; //  which condition are we in
array[N] real phi;
int<lower = 0> L; // number of participants
array[N] int<lower = 1, upper = L> Z;
real prior_theta;
real prior_sigma_u;
real kappa;

}

parameters {
array[K] simplex[8] theta; // mixing proportions - fixed effects

// random effects
array[K, 8] real<lower = 0> sigma_u;
// random effect matrix
array[K, 8, L] real u;

}

transformed parameters {

array[K, 8, L] real u_log_theta;

for (kk in 1:K) {
for (psi in 1:8) {
for (ll in 1:L) {

u_log_theta[kk, psi, ll] = theta[kk][psi] * exp(u[kk, psi, ll]));

}
}
}
}

model {

array[8] real mu;
real piby4 = pi()/4;

mu[1] = 0 * piby4;
mu[2] = 1 * piby4;
mu[3] = 2 * piby4;
mu[4] = 3 * piby4;
mu[5] = 4 * piby4;
mu[6] = 5 * piby4;
mu[7] = 6 * piby4;
mu[8] = 7 * piby4;

// priors

for (kk in 1:K) {

theta[kk] ~ dirichlet(rep_vector(prior_theta, 8));

for (psi in 1:8) {
// we may be better using log_sigma ~ normal
sigma_u[kk, psi] ~ exponential(prior_sigma_u);

for (ll in 1:L) {

u[kk, psi, ll] ~ normal(0, sigma_u[kk, psi]);
}
}
}

int kk, ll; // which condition are we in and which person is it

for (n in 1:N) {

kk = X[n];
ll = Z[n];

vector[8] lps;

for (psi in 1:8) {
// lps[psi] = u_log_theta[kk, psi, ll] + von_mises_lpdf(phi[n] | mu[psi], kappa);
lps[psi] = log(theta[kk, psi]) + von_mises_lpdf(phi[n] | mu[psi], kappa);
}

target += log_sum_exp(lps);
}

}

generated quantities{

// generate priors for easy plotting

simplex[8] pr_theta;
array[8, 5] real<lower = 0> pr_u;
vector[8] pr_sigma_u;

pr_theta = dirichlet_rng(rep_vector(prior_theta, 8));

for (psi in 1:8) {

pr_sigma_u[psi] = exponential_rng(prior_sigma_u);

for (ll in 1:5) {

pr_u[psi, ll] = exp(normal_rng(0, pr_sigma_u[psi]));
}

}
}

``````

Sum of a simplex will always be one; did you mean something else here?

Sorry, I meant theta_u.