Sum to zero constraints and multi-level models - best practise/example code

Modeling
6

Thanks.

I think the missing part of the puzzle is how to integrate this with a variance-covariance matrix.

For example, in my model I have K conditions (usually 2, sometimes 3), and 4 features. Therefore I end up with 4*K parameters (per person).

At present, I define matrix[4*K,L] z_u;:


parameters {

  array[K] real b_a; // weights for class A compared to B  
  array[K] real b_s; // stick-switch rates 
  array[K] real<lower = 0> rho_delta; // distance tuning
  array[K] real rho_psi; // direction tuning

  // random effect variances: 
  // 4*K as we have four fixed effect parameters x K conditions
  vector<lower=0>[4*K] sigma_u;
  cholesky_factor_corr[4*K] L_u;
  // random effect matrix
  matrix[4*K,L] z_u; 
  
}

transformed parameters {

  // this transform random effects so that they have the correlation
  // matrix specified by the correlation matrix above
  matrix[4*K, L] u = diag_pre_multiply(sigma_u, L_u) * z_u;

  // create empty arrays for everything
  array[K] vector[L] u_a, u_stick, u_delta, u_psi;
  // calculate
  for (kk in 1:K) {
    u_a[kk]     = to_vector(b_a[kk]       + u[4*(kk-1)+1]);
    u_stick[kk] = to_vector(b_s[kk]       + u[4*(kk-1)+2]);
    u_delta[kk] = to_vector(rho_delta[kk] + u[4*(kk-1)+3]);
    u_psi[kk]   = to_vector(rho_psi[kk]   + u[4*(kk-1)+4]);
  }
}

I suppose I could try defining array[L] sum_to_zero_vector[4*K] z_u; but I’m not sure how easy it will be to convert z_u to a matrix for the multiplication in transformed parameters {}.

Perhaps the best way forward is to simply remove L_u from the model and skip over estimating all those correlations.

I feel like I am missing something?