How to create factors from factanal?

Question

When performing a factor analysis using factanal the usual result is some loadings table plus several other information. Is there a direct way to use these loadings to create a matrix / data.frame of factors? For example to use them in regression analysis later on.

EDIT: the purpose of this is to obtain variables for subsequent modeling. I only know of factor scores – but suggestions / pointers to other terminology are welcome :)

EDIT2: Joris Meys answer answer is basically what I was asking for. Still though it moves my question towards a direction that might be better suited for statsoverflow, but I will keep it here for now, because the right group of people is the discussing the solution:

What´s the benefit of the regression based scores? The result of the product (ML) is highly correlated with the factors... Honestly I wonder why the difference is that big in my case?

 fa$scores # the correct solution
 fac <- m1 %*% loadings(fa) # the answer on your question
 diag(cor(fac,fa$scores))
 #returns:
Factor1   Factor2   Factor3 
0.8309343 0.8272019 0.8070837 

Answer 1

You asked how to use the loadings for construction of scores. Your solution is, although correct, not doing that. It's using a regression method (alternatively you can use Bartlett's method as well), and this uses the restriction that the scores are uncorrelated, centered around 0 and with variance = 1. These are hence not the same factors as one would obtain by using F = ML with F the factor matrix, M the original matrix and L the loading matrix.

A demonstration with the example from the help files :

v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6)
v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)
v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6)
v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)
v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5)
v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)
m1 <- cbind(v1,v2,v3,v4,v5,v6)

fa <- factanal(m1, factors=3,scores="regression")

fa$scores # the correct solution

fac <- m1 %*% loadings(fa) # the answer on your question

These are clearly different values.

Edit : This has to do with the fact that the Thomson regression scores are based on scaled variables, and take the correlation matrix into account. If you would calculate the scores by hand, you'd do :

> fac2 <- scale(m1) %*% solve(cor(m1)) %*% loadings(fa)
> all.equal(fa$scores,as.matrix(fac2))
[1] TRUE

For more information, see this review

And to show you why it is important : If you calculate the scores the "naive" way, your scores are actually correlated. And that is what you wanted to get rid of in the first place :

> round(cor(fac),2)
        Factor1 Factor2 Factor3
Factor1    1.00    0.79    0.81
Factor2    0.79    1.00    0.82
Factor3    0.81    0.82    1.00

> round(cor(fac2),2)
        Factor1 Factor2 Factor3
Factor1       1       0       0
Factor2       0       1       0
Factor3       0       0       1