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!