Specification of Bayesian SEM models with a data-augmentation approach

Dear Stan folks and @Mauricio_Garnier-Villarre,

I am facing some error messages regardless of the working Stan code for the BSEM (still in the data-augmentation approach).

The stan code I have is:

data{
int N; // sample size
int P; // number of variables
matrix[N, P] Y; // data matrix Y with N rows and P columns
int D; // number of total latent variables
int K; // number of explanatory latent variables (ksi)
int E; // number of outcome latent variables (eta)
matrix[K, K] R; // diagonal matrix
}

parameters{
vector[P] nu; // intercepts of observed variables
vector[P-D] lam; // factor loadings
vector[K] beta; // structural regressions
matrix[N, K] FS_K; // factor scores of explanatory latent variables
vector[E] FS_E; // factor scores of outcome latent variables
vector<lower=0>[P] var_P; // variance for observed variables
cov_matrix[K] cov_K; // covariance for explanatory latent variables
vector<lower=0>[E] var_E; // variance for outcome latent variables
}

transformed parameters{
matrix[N, P] mu_pred; // predicted values of observed variables
vector[N] nu_pred; // predicted values of eta

for(i in 1:N){
mu_pred[i,1] = nu[1] + FS_K[i,1];
mu_pred[i,2] = nu[2] + lam[1]*FS_K[i,1];
mu_pred[i,3] = nu[3] + lam[2]*FS_K[i,1];
mu_pred[i,4] = nu[4] + FS_K[i,2];
mu_pred[i,5] = nu[5] + lam[3]*FS_K[i,2];
mu_pred[i,6] = nu[6] + lam[4]*FS_K[i,2];
mu_pred[i,7] = nu[7] + FS_K[i,3];
mu_pred[i,8] = nu[8] + lam[5]*FS_K[i,3];
mu_pred[i,9] = nu[9] + lam[6]*FS_K[i,3];
mu_pred[i,10] = nu[10] + FS_E[i];
mu_pred[i,11] = nu[11] + lam[7]*FS_E[i];
mu_pred[i,12] = nu[12] + lam[8]*FS_E[i];
}

for(i in 1:N){
  nu_pred[i] = beta[1]*FS_K[i,1]+beta[2]*FS_K[i,2]+beta[3]*FS_K[i,3];
}
}

model{
// hyperpriors
vector[K] u;
u = rep_vector(0, K);

  // priors on intercepts
for (i in 1:P) {nu[i] ~ normal(0, 10);}
  
  // priors on factor loadings
for (i in 1:(P-D)) {lam[i] ~ normal(0, 10);}

// priors on regression coefficients
for (i in 1:K) {beta[i] ~ normal(0, 10);}

// priors on variance
for(i in 1:P) {var_P[i] ~ gamma(1, 0.5);}
var_E ~ gamma(1, 0.5);
cov_K ~ wishart(12, R);

  // likelihood
for(i in 1:N){
  FS_E[i] ~ normal(nu_pred[i], var_E);
  FS_K[i] ~ multi_normal(u, cov_K);
  for(j in 1:P){
    Y[i, j] ~ normal(mu_pred[i,j],var_P[j]);
  }
}
}

The R code I used to run the stan model follows (please note that I saved the Stan code under the name “cfa.stan”:

stan.dat <- list(N=nrow(sim.dat.big), P=ncol(sim.dat.big), K=3, E=1, Y=as.matrix(sim.dat.big),
                 R=matrix(c(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0), nrow=3))

stan.sem.model <- stan_model("cfa.stan")

stan.fit <- rstan::stan(file = "cfa.stan", data = stan.dat, seed = 322, cores = 2, chains = 2,
                        iter = 1000)

What I get as error messages is:

Error in Module(module, mustStart = TRUE) :
Failed to initialize module pointer: Error in FUN(X[[i]], …): no such symbol _rcpp_module_boot_stan_fit4model873ce7b34fc_cfa_mod in package C:/Users/IHN-WH~1/AppData/Local/Temp/RtmpoVLvSt/file873c19b36b87.dll

My questions are:
(1) Do you have any ideas about why these errors are happening (or the meaning of them)?
(2) I did not add the generated quantities{} block for factor indeterminacy (i.e., factor model identification). This is because I already fixed the first factor loadings of each factor to 1 in the Stan code. I am thus wondering whether I still have to do this.
(3) If I have to specify the factor indeterminacy to fix the first factor loadings to 1 instead of fixing the factor variance to 1, how should I specify that? When I referred to the previous Stan issues (Non-convergence of latent variable model) and Stan mailing lists, the method of fixing the factor variance to 1 was used for identification. So, I would like to receive any advice.

Thank you for your time!