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]));
}
}
}