@stijn, thank you again.
Just as you were writing - I was authoring this post!
I did take your suggestion and wrote out the sums in full, decomposing the item family means into an item model mean and a combination of fixed effects:
family_mu = model_mu[mm] + fammu_offset*fam_resid;
family_mu[multten] = family_mu[multten] + char_m10*ident_x;
family_mu[threedig] = family_mu[threedig] + char_d3*ident_y;To accomplish this, as you can see above - I multiplied the scalar values by vectors of 1’s as needed so that I could execute the code without including loops. This also passed the parser - but perhaps my approach (adding ident_x and ident_y ) was unnecessary?
Looking at the code, there is one final modification I would like to make after improving efficiency: I would like to be able to, like the original code, have family-specific random effects. With the modifications here, fam_resid is estimated to be the same across all families. What if I want to estimate vector [M] fam_resid because I think that there is variability across item models in how consistently items were generated? Or if I wanted to estimate vector [P] sigma_item ?
I am currently running three chains - which with the previous version of this same code took… days, with 2000 iterations per chain. It appears already to be running MUCH faster, even with the deep copies. I did already notice in testing versions of this code on simpler models that the non-centered parameterization lead to considerable increases in speed and also I no longer have warnings about low information.