factpca<-function(x=NULL,r=cor(x)){
#Performs Principal Components Factor Analysis using the#
#correlation matrix.  Can be given either the raw data or the#
#correlation matrix.#
xeig<-eigen(r)
xval<-xeig$values
xvec<-xeig$vectors
for (i in 1:length(xval)){
 xvec[,i]<-xvec[,i]*sqrt(xval[i])}
rownames(xvec)<-colnames(r)
return(xvec)
}


factpf<-function(x=NULL,r=cor(x)){
#Performs Principal Factor Factor Analysis using the#
#correlation matrix.  Can be given either the raw data or the#
#correlation matrix.#
n<-ncol(r)
redcormat<-r
diag(redcormat)<-apply(abs(r-diag(1,nrow=n,ncol=n)),2,max)
xeig<-eigen(redcormat)
xval<-xeig$values
xvec<-xeig$vectors
for (i in 1:length(xval[xval>0])){
 xvec[,i]<-xvec[,i]*sqrt(xval[i])}
rownames(xvec)<-colnames(x)
return(xvec[,xval>0])
}


factiter<-function(x=NULL,r=cor(x),niter=100,maxfactors=ncol(r),diag=FALSE){
#Performs Iterated Principal Factor Factor Analysis using the#
#correlation matrix.  Can be given either the raw data or the#
#correlation matrix.#  
n<-ncol(r)
temp<-matrix(0,nrow=n,ncol=n)
comm<-matrix(0,nrow=niter+2,ncol=n)
y<-factpf(r=r)
m<-ncol(y)
temp[1:n,1:m]<-y
comm[1,]<-apply(abs(r-diag(1,nrow=n,ncol=n)),2,max)
comm[2,]<-apply(as.matrix(temp[,1:maxfactors])^2,1,sum)                 
for (i in 3:(niter+2)){
  redcormat<-r
  diag(redcormat)<-comm[i-1,]
  xeig<-eigen(redcormat)
  m<-min(maxfactors,length(xeig$values[xeig$values>0]))
  for (j in 1:m){
     xeig$vectors[,j]<-xeig$vectors[,j]*sqrt(xeig$values[j])} 
  temp<-matrix(0,nrow=n,ncol=n)
  temp[1:n,1:m]<-xeig$vectors[1:n,1:m]
  comm[i,]<-apply(as.matrix(temp[1:n,1:m])^2,1,sum)
}
loadings<-as.matrix(temp[1:n,1:m])
rownames(loadings)<-colnames(r)
rownames(comm)<-(-1:niter)
if (diag!=FALSE)
  {list(comm=comm,loadings=loadings)}
else 
  {loadings}
}



fact<-function(x=NULL,r=cor(x),maxfactors=floor((2*ncol(r)+1-sqrt(8*ncol(r)+5))/2),
	method=c("iter","pf","pca","norm"),rotation=c("none","varimax","quartimax"),
	niter=100,full=FALSE,plot=TRUE,dec=3){
method<-match.arg(method)
rotation<-match.arg(rotation)
if (rotation=="none"){rot="- no rotation"}
if (rotation=="varimax"){rot="- varimax rotation"}
if (rotation=="quartimax"){loaded<-library(GPArotation,logical.return=T)
                           rot<-" - quartimax"
                           if (loaded==F){
                              rot<-"- varimax rotation"
                              rotation<-"varimax"
                              warning("GPArotation library is needed for",
                                " quartimax rotation.  Varimax was used instead",
                                 call.=FALSE)}}
if (method=="iter"){ftemp<-factiter(r=r,niter=niter,maxfactors=maxfactors,diag=T)
		    loadings<-ftemp$loadings
                    concheck<-max(abs(ftemp$comm[niter+2,]-ftemp$comm[niter+1,]))
		    if (concheck>=5e-4){
                         warning("Convergence not achieved, difference of ",
                         as.character(round(concheck,5)),
			 " after ",as.character(niter)," iterations.",
                         call.=FALSE)}
		    maxcon<-apply(ftemp$comm,1,max)
                    maxlegal<-max(as.integer(names(maxcon[maxcon<1])))
                    if (as.integer(maxlegal)<niter){
                         warning("Communality greater than one, ",maxlegal,
                         " was last legal iteration.",call.=FALSE)}
                    meth<-paste("iterated principal factor",rot)}
if (method=="pf"){loadings<-factpf(r=r)
                 loadings<-loadings[,1:min(ncol(loadings),maxfactors)]
		  meth<-paste("principal factor",rot)}
if (method=="pca"){loadings<-factpca(r=r)
                  loadings<-loadings[,1:min(ncol(loadings),maxfactors)]
		  meth<-paste("principal components",rot)}
if (method=="norm"){q<-ncol(r)
		   maxfactors<-min(maxfactors,
                         floor((2*ncol(r)+1-sqrt(8*ncol(r)+5))/2))
		   faout<-factanal(covmat=r,factors=maxfactors,rotation="none")
                   loadings<-diag(rep(1,ncol(r)))%*%faout$loadings
                   meth<-paste("multivariate normal",rot)}
loadings<-as.matrix(loadings)
if ((rotation=="varimax")&&(ncol(loadings)>1)){
   loadings<-diag(rep(1,ncol(r)))%*%varimax(loadings)$loadings}
if ((rotation=="quartimax")&&(ncol(loadings)>1)){
   loadings<-quartimax(loadings,normalize=TRUE)$loadings}
loadings<-as.matrix(loadings)
loadings<-loadings[,order(apply(-loadings^2,2,sum))]
loadings<-as.matrix(loadings)
colnames(loadings)<-paste("Factor",as.character(1:ncol(loadings)),sep="")
rownames(loadings)<-colnames(r)
resid<-r-loadings%*%t(loadings)-diag(1-apply(loadings^2,1,sum))
cortri<-r[upper.tri(r)]
restri<-resid[upper.tri(resid)]
predtri<-cortri-restri
if (plot==TRUE){
  par(mfrow=c(2,2))
  plot(princomp(covmat=r),main="")
  plot(cortri,predtri,xlab="Original Correlation",
	ylab="Predicted Correalation",xlim=c(-1,1),ylim=c(-1,1))
  lines(c(-1,1),c(-1,1))
  plot(predtri,restri,xlab="Predicted Correlation",ylab="Resdisual Correlation",
	main=paste("r-squared=",as.character(round(cor(predtri,restri)^2,3))))
  lines(c(-1,1),c(0,0))
  abline(lm(restri~predtri))
  hist(restri,xlab="Residual Correlations",
        main=paste("mean=",as.character(round(mean(restri),3)),
                   "  s.d.=",as.character(round(sd(restri),3))))
  par(mfrow=c(1,1))
}
loadings<-t(t(loadings)*sign(apply(sign(loadings),2,sum)+.001))
evs<-eigen(r)$values
comm<-apply(loadings^2,1,sum)
names(comm)<-colnames(r)
vexpl<-apply(loadings^2,2,sum)
vexpl<-rbind(vexpl,vexpl/ncol(r))
colnames(vexpl)<-colnames(loadings)
rownames(vexpl)<-c("variance explained","percent explained")
resid<-resid
if (full!=TRUE){resid<-summary(resid[upper.tri(resid)])}
list(eigen.values=round(evs,dec),
      method=meth,
      loadings=round(loadings,dec),
      communalities=round(comm,dec),
      importance=round(vexpl,dec),
      residuals=round(resid,dec))
}