Mixed effects von Mises models using a unit vector?

TLDR

  1. I have working stan code, but would like to be sure that it is actually modelling: mu ~ 1 + (1 | animal) on a circle, as well as log(kappa) ~ 1 + (1 | animal).
  2. I would like to know how to make a more general family of von Mises models in brms where mu is estimated as a unit_vector, using custom_family, stanvars etc.

Details
I’ve now had some success with implementing the unit_vector transform as a link function for mixed-effects models using the von Mises distribution.


Unfortunately, I just haven’t been able to get my head around post processing the Stan model, or the appropriate syntax for adding this transform to a brms custom_family. I’ve tried reverse engineering brms generated code handwritten.uvMEmises1.stan (3.2 KB) in the hope that I can work out how to define the custom_family. The main problem seems to be that if I follow the original post faithfully, then this requires the introduction of a new parameter to be estimated within the main modelling loop, not just a new link function and likelihood function. Any advice would be much appreciated!

In the original post the lpdf was calculated for mu_vec as well as raw angles (should be the same I guess?), so I’ve added this into the optimisation loop:

  // likelihood including all constants
    if (!prior_only) {
      for (n in 1:N) {
        target += von_mises_lpdf(Y[n] | mu[n], kappa[n]);
        //add lpdf from the unit_vector estimate
        //I am not sure how Stan updates mu_vec
        target += von_mises_unitvector_lpdf(mu_vec | mu[n], kappa[n]);
      }
    }

Where von_mises_unitvector_lpdf is the unit_vector version of the default brms von_mises_real_lpdf (which I was able to insert as a stavar).
Using custom_families I can only insert this into the modelling block as a separate loop outside of the main modelling loop unit vector MEvonMises brms attempt.R (7.7 KB) , since stanvars inserts code between the parameter declaration and modelling block by default. This is particularly problematic when kappa is estimated using the log link (which is generally a good idea) since this link is applied below the unit_vector loop I can insert.

It would be great to have a brms method for modelling this for any number of animals and fixed effects. My colleagues and I work with this kind of dataset a lot and I would be a lot more comfortable using brms for fitting and post-processing.