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:
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?
Are the results from one of the two more valid than the other, i.e. either Stata or R psych?
Is there a recommendation to use principal components, pca
psych::principal
or principal component factorspcf
-- for a polychoric correlation matrix in R psych, lavaan, etc.? (Stata-folks seem to agree onpcamat
~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.