I need to build a CFA model that will be then used as a component of a more complex model. Since the final model fairly big, I need the CFA to be as fast and efficient as possible.
I built two CFA models: one with data augmentation, and another that draws from the marginal (something like what blavaan does). Although the blavaan-type model should run faster, it’s like 2.5 times slower. Is there any way I could speed it up?
Thanks a lot for your help!
The marginal model:
data {
int<lower=1> N; // n of individuals/observations
int<lower=1> M; // n of factors/latent variables
int<lower=1> P; // n of observed indicators
vector[P] Ind[N]; // matrix[N,P] of indicators
int<lower=0, upper=1> D_skeleton[P,M]; // which LVs are related to which indicators? 1 if yes, 0 if not
int Dvec_length; // equal to sum(D_skeleton)
}
parameters {
vector<lower=0>[Dvec_length] Dvec; // factor loadings
vector<lower=0>[P] psi; // variance of indicators
vector[P] d; // mean of indicators
cholesky_factor_corr[M] L_Omega; // cholesky factor or the covariance of LVs. (Actually a correlation so the variance of the LVs = 1)
}
transformed parameters {
matrix[P,M] D = Dvec_to_D(P, M, D_skeleton, Dvec); // convert to a matrix of factor loadings. See the function below.
corr_matrix[M] Omega = multiply_lower_tri_self_transpose(L_Omega);
}
model {
matrix[P,P] new_cov = quad_form(Omega, D') + diag_matrix(psi);
Dvec ~ normal(1, 2);
psi ~ inv_gamma(1,2);
d ~ normal(0, 10);
L_Omega ~ lkj_corr_cholesky(1);
Ind ~ multi_normal(d, new_cov);
}
The mathematical model:
I_i = \delta + D \nu_i + \eta_i \quad \forall i
\nu_i \sim N_M(0, \Omega) (the latent variables)
\eta_i \sim N_P(0, \Psi) (the measurement errors)
So the reduced form model is:
I_i \sim N_P(\delta, D\Omega D^T + \Psi)
The Dvec_to_D
function:
functions {
// Transform a vector of factor loadings into a matrix with appropriate zero constraints
// @param P the number of observed indicators
// @param M the number of factors/latent variables
// @param D_skeleton a matrix with 1s where factor loadings should be estimated, 0s otherwise
// @param Dvec a vector with factor loadings
// @returns a [P,M] matrix with zeroes where D_skeleton is equal to zero
matrix Dvec_to_D(int P,
int M,
int[ , ] D_skeleton,
vector Dvec) {
matrix[P,M] ret_matrix;
int counter = 1;
for (i in 1:P) {
for (j in 1:M) {
if (D_skeleton[i,j] == 1) {
ret_matrix[i,j] = Dvec[counter];
counter += 1;
}
else {
ret_matrix[i,j] = 0;
}
}
}
return ret_matrix;
}
}