[R] Factor Analysis by Principal Components (Code Attached)
Brett Magill
bmagill at earthlink.net
Thu Aug 1 19:48:44 CEST 2002
Hi all,
I have written some (very rudimentary) code to produce "factor analysis"
results using the principal components method.
I would appreciate it if anyone who is interested in this could could try it
out and compare it to results with other packages given known results, which I
have little ability to do. If other find it useful, I could clean it up, make
it more general, add functionality, etc.
With respect to interpretation, I have found my results in line with
factanal(), at least with respect to interpretation, but generally higher. I
compared a couple of examples with some SAS results where I could find data
and results on the web--again, comparable. There seems to be some difference
between the principal components estimated by the two packages???
I use princomp to get the loadings and the eigen values. The fa coefficients
are then the PC loadings * the sqrt of the eigenvalues (from Johnson and
Wichern) for each retained factor. These coeffcients are then rotated with
varimax().
The function takes a covariance or correlation matrix as its input. If you do
not supply an n for factors, it uses the eigen > 1 criterion.
-----------Code below-----------------
function (xmat, factors=NULL, cor=TRUE) {
prc <- princomp ( covmat = xmat ,cor = cor )
eig <- prc$sdev^2
if (is.null(factors)) factors <- sum ( eig >= 1 )
loadings <- prc$loadings [ , 1:factors ]
coefficients <- loadings [ , 1:factors ] %*% diag ( prc$sdev[1:factors] )
rotated <- varimax ( coefficients ) $ loadings
fct.ss <- apply( rotated, 2 , function (x) sum (x^2) )
pct.ss <- fct.ss / sum (eig)
cum.ss <- cumsum ( pct.ss )
ss <- t ( cbind ( fct.ss , pct.ss, cum.ss ) )
return ( coefficients , rotated , ss )
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list