I agree that the sufficient statistic approach can be faster, especially for large N.
Another thing to consider is marginalizing over eta. Instead of Y[, k] ~ normal()
, you would have something like
psi = quad_form_sym(Omega, diag_matrix(tau));
for (i in 1:N) {
Y[i, ] ~ multi_normal(rep_vector(0, K), quad_form_sym(psi, lambda'));
}
You might also look at the blavaan R package, which contains many of these optimizations and uses lavaan syntax for model specification. Here is an old post linking to some other related information: