Agree, the blavaan exported Stan code, is bulky and can be difficult to follow. I still like to recommend these options, as its a good way to learn coding tricks by seeing code that works
As of the model, I worked and example model similar to the original model posted, with the corrected signed parameters in the generated quantities block. Some comments about the model, I estimated the model based on the data augmentation approach, which for me at least make it easier to extend for Multilevel CFA. But the log-likelihood for model comparison is estimated with the covariance matrix approach, so the marginal log-likelihood is used instead of the conditional, as recommended by
Merkle, E. C., Furr, D., & Rabe-Hesketh, S. (2019). Bayesian model assessment: Use of conditional vs marginal likelihoods. Psychometrika , 84 , 802–829.
data{
int N; // sample size
int P; // number of variables
int D; // number of dimensions
vector[P] X[N]; // data matrix of order [N,P]
int n_lam; // how many factor loadings are estimated
}
parameters{
vector[P] b; // intercepts
matrix[N,D] FS_UNC; // factor scores, matrix of order [N,D]
cholesky_factor_corr[D] L_corr_d; // cholesky correlation between factors
vector<lower=0>[P] sd_p; // residual sd for each variable
vector<lower=-30, upper=30>[n_lam] lam_UNC; // not positivily constraint
}
transformed parameters{
vector[D] M; // factor means
vector<lower=0>[D] Sd_d; // sd for each factor
vector[P] mu_UNC[N]; // predicted values
cholesky_factor_cov[D] L_Sigma;
M = rep_vector(0, D);
Sd_d = rep_vector(1, D);
L_Sigma = diag_pre_multiply(Sd_d, L_corr_d);
for(i in 1:N){
mu_UNC[i,1] = b[1] + lam_UNC[1]*FS_UNC[i,1];
mu_UNC[i,2] = b[2] + lam_UNC[2]*FS_UNC[i,1];
mu_UNC[i,3] = b[3] + lam_UNC[3]*FS_UNC[i,1];
mu_UNC[i,4] = b[4] + lam_UNC[4]*FS_UNC[i,2];
mu_UNC[i,5] = b[5] + lam_UNC[5]*FS_UNC[i,2];
mu_UNC[i,6] = b[6] + lam_UNC[6]*FS_UNC[i,2];
mu_UNC[i,7] = b[7] + lam_UNC[7]*FS_UNC[i,3];
mu_UNC[i,8] = b[8] + lam_UNC[8]*FS_UNC[i,3];
mu_UNC[i,9] = b[9] + lam_UNC[9]*FS_UNC[i,3];
}
}
model{
L_corr_d ~ lkj_corr_cholesky(1);
b ~ normal(0, 10);
lam_UNC ~ normal(0, 10);
sd_p ~ cauchy(0,2.5);
///Sd_d ~ cauchy(0,2.5);
for(i in 1:N){
for(j in 1:P){
X[i,j] ~ normal(mu_UNC[i,j], sd_p[j]);
}
FS_UNC[i] ~ multi_normal_cholesky(M, L_Sigma);
}
}
generated quantities{
real log_lik[N]; ///// log likelihood matrix
real dev; /////// global deviance
real log_lik0; ///// global log likelihood
vector[N] log_lik_row;
cov_matrix[P] Sigma;
cov_matrix[D] Phi_lv;
matrix[P,P] lambda_phi_lambda;
cholesky_factor_cov[P] L_Sigma_model;
matrix[P,P] theta_del;
corr_matrix[D] Rho_UNC; /// correlation matrix
corr_matrix[D] Rho; /// correlation matrix
matrix[P,D] lam; // factor loadings
matrix[N,D] FS; // factor scores, matrix of order [N,D]
///vector[P] mu[N]; // predicted values
/// signed corrected parameters
Rho_UNC = multiply_lower_tri_self_transpose(L_corr_d);
Rho = Rho_UNC;
FS = FS_UNC;
lam = rep_matrix(0, P, D);
lam[1:3,1] = to_vector(lam_UNC[1:3]);
lam[4:6,2] = to_vector(lam_UNC[4:6]);
lam[7:9,3] = to_vector(lam_UNC[7:9]);
// factor 1
if(lam_UNC[1] < 0){
lam[1:3,1] = to_vector(-1*lam_UNC[1:3]);
FS[,1] = to_vector(-1*FS_UNC[,1]);
if(lam_UNC[4] > 0){
Rho[1,2] = -1*Rho_UNC[1,2];
Rho[2,1] = -1*Rho_UNC[2,1];
}
if(lam_UNC[7] > 0){
Rho[1,3] = -1*Rho_UNC[1,3];
Rho[3,1] = -1*Rho_UNC[3,1];
}
}
// factor 2
if(lam_UNC[4] < 0){
lam[4:6,2] = to_vector(-1*lam_UNC[4:6]);
FS[,2] = to_vector(-1*FS_UNC[,2]);
if(lam_UNC[1] > 0){
Rho[2,1] = -1*Rho_UNC[2,1];
Rho[1,2] = -1*Rho_UNC[1,2];
}
if(lam_UNC[7] > 0){
Rho[2,3] = -1*Rho_UNC[2,3];
Rho[3,2] = -1*Rho_UNC[3,2];
}
}
// factor 3
if(lam_UNC[7] < 0){
lam[7:9,3] = to_vector(-1*lam_UNC[7:9]);
FS[,3] = to_vector(-1*FS_UNC[,3]);
if(lam_UNC[1] > 0){
Rho[3,1] = -1*Rho_UNC[3,1];
Rho[1,3] = -1*Rho_UNC[1,3];
}
if(lam_UNC[4] > 0){
Rho[3,2] = -1*Rho_UNC[3,2];
Rho[2,3] = -1*Rho_UNC[2,3];
}
}
/// marginal log-likelihood based on signed corrected parameters
Phi_lv = quad_form_diag(Rho, Sd_d);
lambda_phi_lambda = quad_form_sym(Phi_lv, transpose(lam));
theta_del = diag_matrix(sd_p);
Sigma = lambda_phi_lambda + theta_del;
L_Sigma_model = cholesky_decompose(Sigma);
for(i in 1:N){
log_lik[i] = multi_normal_cholesky_lpdf(X[i] | b, L_Sigma_model);
}
log_lik0 = sum(log_lik); // global log-likelihood
dev = -2*log_lik0; // model deviance
}
I compared this model with lavaan and blavaan, with simulated data
library(lavaan)
library(loo)
library(blavaan)
library(rstan)
#rstan_options(auto_write = TRUE)
#options(mc.cores = parallel::detectCores())
mod1 <-'
f1 =~ .5?y1+.5?y2+.5?y3 #+.1?y4 +.2?y8
f2 =~ .4?y4+.4?y5+.4?y6 #+.1?y1 +.2?y7
f3 =~ .6?y7+.6?y8+.6?y9 #+.1?y6 +.2?y3
f1 ~~ .3?f2+.4?f3
f2 ~~ .5?f3'
mod2 <-'
f1 =~ y1+y2+y3
f2 =~ y4+y5+y6
f3 =~ y7+y8+y9
f1 ~~ f2+f3
f2 ~~ f3'
datsim <- simulateData(mod1,sample.nobs=300,std.lv=T)
head(datsim)
fit <- cfa(mod2,data=datsim,std.lv=T,meanstructure=T)
summary(fit,standardized=T,fit.measures=T)
fitb <- bcfa(mod2,data=datsim,std.lv=T,target="stan")
summary(fitb,standardized=T)
fitMeasures(fitb)
###### run the Stan model
X <- as.matrix(datsim)
P <- ncol(datsim)
N <- nrow(datsim)
D <- 3
dat <- list(N=N,P=P,D=D,X=X, n_lam=9)
param <- c("b","lam","Rho","sd_p","M","Sd_d",
"log_lik0","dev","log_lik")#
##### loop to check for convergence
ta2 <- Sys.time()
iters <- 5000
keep <- 3000
rhat <- 20
while(rhat > 1.05){
iters <- iters + 3000
burn <- iters - keep
OUT_st <- stan(data=dat, pars=param,
model_code=cfa_stan_fv,
chains=3, iter=iters,warmup=burn,
control=list(adapt_delta=.99, max_treedepth=20))
ss <- summary(OUT_st)$summary
print(rhat <- max(ss[,"Rhat"], na.rm=T))
}
tb2 <- Sys.time()
tb2 - ta2
OUT_st
print(OUT_st,digits=3, pars=c("b","lam","Rho","sd_p","M","Sd_d",
"log_lik0","dev"))
log_lik1 <- extract_log_lik(OUT_st, merge_chains = FALSE)
(waic_cfa <- waic(log_lik1))
(loo_cfa <- loo(log_lik1, r_eff=relative_eff(exp(log_lik1))))
I teach a BSEM summer course, havent worked a Multilevel model yet, I will seek to develop a couple of useful examples from this. Its an interesting issue
