Is there a way to obtain factor correlations with significance levels in R's lavaan for CFA?

I’m working on a Confirmatory Factor Analysis (CFA) for my survey using R’s lavaan package. I’ve managed to get the correlation matrix for my factors, but I’m stuck on how to include the p-values for each correlation. I’ve tried using cov2cor(inspect(fit, what = "est")$psi), but it only gives me the matrix without significance levels.

Here’s a basic example of my CFA model:

my_cfa <- '
factor1 =~ q1 + q2 + q3
factor2 =~ q4 + q5 + q6
factor3 =~ q7 + q8 + q9
'

results <- cfa(my_cfa, data = my_dataset,
               auto.var = TRUE, auto.fix.first = TRUE,
               auto.cov.lv.x = TRUE)

Any tips on how to get both the correlations and their p-values in one output? It would really help with interpreting my results. Thanks!

Hey there, SurfingWave! :wave:

I totally get your struggle with lavaan - it’s a powerful package, but sometimes getting exactly what you want can be tricky!

Have you tried using the parameterEstimates() function? It might give you what you’re looking for. Something like this could work:

pe <- parameterEstimates(results, standardized = TRUE)
factor_correlations <- pe[pe$op == "~~" & pe$lhs != pe$rhs, ]

This should give you a nice table with the factor correlations, standard errors, and p-values.

But I’m curious - what’s the end goal of your analysis? Are you looking to validate a specific theoretical model, or is this more exploratory? And how many factors are you working with in total?

It’d be cool to hear more about your research. Maybe there are some other approaches we could explore depending on what you’re trying to achieve!

yo surfingwave, i’ve had this prob too. try the lavTestLRT() function. it gives likelihood ratio tests for factor correlations. like this:

full_model <- cfa(my_cfa, data=my_dataset)
reduced_model <- cfa(my_cfa, data=my_dataset, orthogonal=TRUE)
lavTestLRT(full_model, reduced_model)

this compares ur model to one where factors arent correlated. p-values show if correlations r significant. hope this helps!

I’ve encountered a similar issue in my research. One approach that worked for me was using the modindices() function in combination with parameterEstimates(). Here’s what I did:

pe <- parameterEstimates(results, standardized = TRUE)
factor_cors <- pe[pe$op == '~~' & pe$lhs != pe$rhs, ]
mi <- modindices(results)
factor_cors_with_p <- merge(factor_cors, mi, by.x = c('lhs', 'rhs'), by.y = c('lhs', 'rhs'))

This gives you a dataframe with correlations, standard errors, and p-values. It’s not the most elegant solution, but it gets the job done. Just be cautious when interpreting these p-values, especially with larger sample sizes, as they can be overly sensitive.

Have you considered using a Bayesian approach? It might offer a more nuanced interpretation of factor relationships, especially if you’re dealing with complex models or limited sample sizes.