rock <- read.csv("http://www.stat.ufl.edu/~winner/data/rockstrength.csv") attach(rock); names(rock) X <- as.matrix(rock[,2:10]) S <- cov(X) R <- cor(X) n <- nrow(X) rock.fa <- lapply(1:4, function(nf) factanal(rock[,2:10],factors=nf,method="mle",scores="regression")) rock.fa (scores <- rock.fa[[4]]$scores) ## Difference between R and LL'+psi pred <- rock.fa[[4]]$loadings %*% t(rock.fa[[4]]$loadings) + diag(rock.fa[[4]]$uniquenesses) R-pred par(pty="s") plot(rock.fa[[4]]$scores[,1],rock.fa[[4]]$scores[,2], ylim=range(rock.fa[[4]]$scores[,1]),xlab="Factor1", ylab="Factor2")