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.