p = 100
data = matrix(scan("BCS-males.txt"),byrow=TRUE,ncol=p)
y1 = data[,c(1:7)]
y2 = data[,47+c(2,17,19,20,23,30,31,33,52,53)]
y  = cbind(y1,y2)
names1 = c("Picture Language","Friendly Math","Reading","Matrices","Recall digits",
           "Similarities","Word definition")
names2 = c("Child Dev 2","Child Dev 17","Child Dev 19","Child Dev 20","Child Dev 23",
           "Child Dev 30","Child Dev 31","Child Dev 33","Child Dev 52","Child Dev 53")
names = c(names1,names2)           
p = ncol(y)
n = nrow(y)
round(cor(y))

y = y - matrix(apply(y,2,mean),n,p,byrow=TRUE)
y = y%*%diag(1/sqrt(apply(y,2,var)))

pdf(file="histograms.pdf",width=10,height=8)
par(mfrow=c(4,4))
for (i in 1:16)
  hist(y[,i],xlab="",prob=TRUE,main=names[i])
dev.off()

pdf(file="pairs1.pdf",width=10,height=8)
pairs(y1,label=names1,pch=16)
dev.off()
pdf(file="pairs2.pdf",width=10,height=8)
pairs(y2[,1:7],label=names2[1:7],pch=16)
dev.off()


fac0 = factanal(y,2,rotation="none",score="regression")
cbind(names,round(cbind(-fac0$load[,1],fac0$load[,2],fac0$uniq),2))

fac1 = factanal(y,2,rotation="varimax",score="regression")
cbind(names,round(cbind(-fac1$load[,1],fac1$load[,2],fac0$uniq),2))

fac1 = factanal(y,2,rotation="none",score="regression")
fac2 = factanal(y,2,rotation="varimax",score="regression")

pdf(file="varimax.pdf",width=10,height=8)
par(mfrow=c(1,1))
plot(fac1$load,xlim=range(fac1$load[,1],fac2$load[,1]),ylim=range(fac1$load[,2],fac2$load[,2]),pch=16,cex=2,
     xlab="First common factor",ylab="Second common factor")
abline(h=0);abline(v=0)
for (i in 1:p) segments(fac1$load[i,1],fac1$load[i,2],fac2$load[i,1],fac2$load[i,2],col=grey(0.5),lwd=3)
for (i in 1:p) points(fac1$load[i,1],fac1$load[i,2],pch=16,cex=2)
for (i in 1:p) points(fac2$load[i,1],fac2$load[i,2],col=2,pch=16,cex=2)
legend(0.4,0.8,legend=c("FA","FA+VARIMAX"),col=1:2,pch=16,bty="n")
dev.off()