Fitting of a Bayes CFA Model in stan

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:

1 Like