Howdy howdy, I had a quick look at this.
spp_intercept ~ normal( 0 , sppvar );
The second argument to normal for a Stan model is a standard deviation, but that’d only mess up your labeling.
matrix[N_site_id, D] latent_vars;
These are the zs, right? Aren’t there just N_sites of them?
The matrix uij is a function of z and lambda. Given those are both parameters, it’s probably going to be painful to sample that (from the paper “Estimation is difficult because neither the latent variables zi nor the factor loadings lj are known”). How did they get around this?
It’s probably worth your time to generate stuff from the approximate model exactly. This’ll help you isolate issues with the model (if you know exactly how data was generated – hopefully you can get it to fit). This latent variable vs. factor loadings thing seems like something worth understanding.
If that works out, then you can move to working with this model when it’s an approximation.