I have some trouble translating between Stata and R (psych, lavaan) about factors and components with an interest in the loadings/eigenvectors of the items. I seek to validate results (understand the results I am getting, and getting similar results) and ultimately choose the correct method in the software. The terminology is somewhat difficult. And other SE users have found it difficult to compare results with Stata.

In Stata, it is possible to feed polychoric correlations to factormat to "perform a factor analysis... principal factor, iterated principal factor, principal-component factor, and maximum-likelihood factor analyses". Principal-component factor, pcf--similar to pca except using different normalizations; to unit length or eigenvalue, as also noted on SE--is recommended by Acock: "when you're trying to develop a measure of a concept". Stas Kolenikov uses something similar in his polychoricpca. Indeed, thanks to StasK polychoric is possible in Stata. rotate "performs a rotation of the loading matrix". An example:

polychoric item1 item2 item3..., pw // pairwise
matrix r = r(R)
local N = r(sum_w)
factormat r, n(`N') pcf // similar to    _pcamat r, n(`N')_
rotate, promax

Results I get make sense. And they compare with results in the sem framework, also when using adf estimates, and ologit functions.

Turning to R, I am not the only who has experienced issues feeding polychoric correlation matrix to factor analyses using psych, even when specifying other scoring methods than "regression" as ttnphns seems to suggest. ttnphns: "I don't know psych and its options but I suspect that the package just will process such matrix as if it were Pearson correlation." psych::factor.scores:

factor.scores uses four different ways of estimate factor scores. In all cases, the factor score estimates are based upon the data matrix, X, times a weighting matrix, W, which weights the observed variables.

For polytomous or dichotomous data, factor scores can be estimated using Item Response Theory techniques (e.g., using link{irt.fa} and then link{scoreIrt}. Such scores are still just factor score estimates, for the IRT model is a latent variable model equivalent to factor analysis.

# sparing you for detailed variants of arguments to fa()

efa <- psych::fa(data, cor = "poly", 
          scores = "Bartlett",
          rotate = "promax")
# output
In cor.smooth(R) : Matrix was not positive definite, smoothing was done

efa$r # shows the correlation matrix is quite off or misguided.

psych::polychoric(data) # gives similar warning and result

Results make a lot less sense. pca psych::principal showed similar less meaningful results. With Pearson's correlation matrix it all makes much more sense.

  • It does warn of eigenvector (principal components) smoothing where "Negative eigen values are replaced with 100 * eig.tol, the matrix is reproduced and forced to a correlation matrix using cov2cor."
  • The correlation matrix from polychoric appear to be incorrect.

I guess my question is three-fold:

  1. What are, or could be, the reasons for diverging results from procedures to feed polychoric correlation matrices to factor/components analyses in Stata and R psych provided above? - only the scoring procedure?

  2. Are the results from one of the two more valid than the other, i.e. either Stata or R psych?

  3. Is there a recommendation to use principal components, pca psych::principal or principal component factors pcf -- for a polychoric correlation matrix in R psych, lavaan, etc.? (Stata-folks seem to agree on pcamat~factormat pcf.)

Note: General, non-software-specific problems with feeding polychoric matrices to f/c analyses

According to ttnphns, "when the analysis is based on polychoric correlations, the traditional methods of computation of factor scores are inapplicable". ttnphns has also explained:

One disadvantage of the approach ["PCA/FA performed on polychoric (for ordinal data) correlations"] is that - at inferring the correlations - it has no clues to the multivariate distribution of the underlying variables, - can "conceive of" at most bivariate distributions, thus bases itself not on full information.

The argument is tightened in the comments below:

The problem is not with extraction of loadings from the correlations, but with the data to use, to compute the scores. The correlations and hence the loadings do not correspond to the data directly anymore.

  • 1
    $\begingroup$ Unfortunately, I (whom you've cited) cannot help you with a recipe. I can confirm and repeat that classic methods used to compute f. scores in linear FA would be inadequate in the case of tetra/polychoric correlations input. The problem is not with extraction of loadings from the correlations, but with the data to use, to compute the scores. The correlations and hence the loadings do not correspond to the data directly anymore. I know special methods exist, but I've never used them and don't know a package providing them. Lavaan? Psych? Maybe, but I'm not an R expert at all. $\endgroup$
    – ttnphns
    Commented Jun 8, 2023 at 21:56
  • 1
    $\begingroup$ Once again: the problem is not with factor extraction stage. So, you doubt about using pca or pcf is irrelevant to the problem, for me. $\endgroup$
    – ttnphns
    Commented Jun 8, 2023 at 22:01
  • $\begingroup$ Thanks so much @ttnphns, I will rephrase my question. But do you have any clue why using this technique could differ in Stata vs R/psych, knowing the translation issues that have already been reported? $\endgroup$
    – Johan
    Commented Jun 9, 2023 at 5:27
  • 2
    $\begingroup$ Check stepwise, onе аt a time, comparing Stata and R results. Are computed correlation matrices the same? Are extracted by the same method loadings the same? Are rotated by the same method loadings the same? Are f. scores computed by the same method the same? $\endgroup$
    – ttnphns
    Commented Jun 9, 2023 at 10:45
  • $\begingroup$ More generally, as algorithms diverge from linear to more nonlinear routines, the interpretation of the results becomes increasingly less intuitive, making less sense. This is my experience. $\endgroup$
    – user78229
    Commented Jun 9, 2023 at 10:58

1 Answer 1


It turned out the correlation matrix of psych::fa() and psych::polychoric() was off because I had rescaled the items to 0-1, and psych apparently assumed otherwise without raising an error. So the correlations were weird and probably incorrect.

The POLYCHORIC_R() function in the R package EFA.dimensions produced a matrix just fine. At least it looked more similar to what everything else produces. It initially tried to use psych::polychoric() by Revelle.

The data matrix does not appear to consist of whole numbers and is therefore not appropriate for the computation of polychoric correlations.

Consider stopping the program.

A warning occurred when using method=="Revelle". method was therefore changed to "Fox"

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.