How to get the same results with factanal and psych::fa in factor analysis with R
There are two common options to run a factor analysis in R, either with the base factanal
function or the fa
function from the package psych
. The base function factanal
uses maximum likelihood estimation, while the fa
function from the psych
package uses minimum residual estimation as the default. Is it possible to reproduce a factanal result with the fa function? Yes, one just has to set the correct parameters.
As an example let us fit the big5 model. First, we prepare the data and libraries:
library(psych)
library(qgraph) # for big5 data set
data(big5)
Now, we can fit the two models. factanal
uses varimax rotation as default, so we also need to set this for fa
.
fa1 <- factanal(big5, 5)
fa2 <- fa(big5, 5, fm = "ml", rotate = "varimax")
loadings1 <- unclass(fa1$loadings)
loadings2 <- unclass(fa2$loadings)
Now let us look at the difference between the factor loadings. We round the results to three decimal points:
diff <- round(loadings1, 3) - round(loadings2, 3)
mean(diff > 0)
## [1] 0
There is no difference, so fa
can be used instead of factanal
. fa
is much more flexible, so I would recommend to always use it. For instance, it offers plenty of rotation options:
orthogonal:
“none”, “varimax”, “quartimax”, “bentlerT”, “equamax”, “varimin”, “geominT” and “bifactor” are orthogonal rotations.
oblique:
“Promax”, “promax”, “oblimin”, “simplimax”, “bentlerQ,”geominQ” and “biquartimin” and “cluster”
There are also many factoring options:
“minres”, “uls”, “ols”, “wls”, “gls”, “pa”, “ml”, “minchi”, “minrank”, “alpha”
With the function fa.parallel
we get a parallel analysis with a plot:
options(mc.cores = 1)
pa <- fa.parallel(big5, nfactors = 5, fm = "ml")
## Loading required namespace: GPArotation
## Parallel analysis suggests that the number of factors = 21 and the number of components = 16
To better see where the empirical eigen values are larger than the 95% quantile of the simulated ones, we can also make a table:
head(data.frame(pa$fa.values, pa$fa.sim), 30)
## pa.fa.values pa.fa.sim
## 1 19.9050107 1.8485375
## 2 11.2511343 1.7748310
## 3 10.9448784 1.7266967
## 4 8.0589502 1.6849484
## 5 5.5699333 1.6428898
## 6 4.5693196 1.5888678
## 7 3.4607269 1.5486704
## 8 3.0821743 1.5198653
## 9 2.6852770 1.4891862
## 10 2.5276525 1.4572315
## 11 2.2909136 1.4261564
## 12 1.9957009 1.3962894
## 13 1.7590163 1.3705031
## 14 1.6678141 1.3408697
## 15 1.5389754 1.3139674
## 16 1.4745543 1.2889220
## 17 1.4650375 1.2647211
## 18 1.4350744 1.2457785
## 19 1.3307514 1.2209643
## 20 1.2633818 1.1933441
## 21 1.2150480 1.1754355
## 22 1.1498082 1.1516539
## 23 1.1272836 1.1339377
## 24 1.0908419 1.1118007
## 25 1.0703017 1.0889166
## 26 1.0205069 1.0703950
## 27 1.0023410 1.0532331
## 28 0.9782625 1.0324438
## 29 0.9585486 1.0096161
## 30 0.9254454 0.9911679
Quite handy. In conclusion, use psych::fa
and psych::fa.parallel
instead of factanal
!