"principal" <- function(r,nfactors=1,residuals=FALSE,rotate="varimax",n.obs = NA, scores=FALSE,missing=FALSE,impute="median") { cl <- match.call() n <- dim(r)[2] if (n!=dim(r)[1]) { n.obs <- dim(r)[1] if(scores) {x.matrix <- r if(missing) { miss <- which(is.na(x.matrix),arr.ind=TRUE) if(impute=="mean") { item.means <- colMeans(x.matrix,na.rm=TRUE) #replace missing values with means x.matrix[miss]<- item.means[miss[,2]]} else { item.med <- apply(x.matrix,2,median,na.rm=TRUE) #replace missing with medians x.matrix[miss]<- item.med[miss[,2]]} }} r <- cor(r,use="pairwise") # if given a rectangular matrix, then find the correlations first } else { if(!is.matrix(r)) { r <- as.matrix(r)} sds <- sqrt(diag(r)) #convert covariance matrices to correlation matrices r <- r/(sds %o% sds) } #added June 9, 2008 if (!residuals) { result <- list(values=c(rep(0,n)),rotation=rotate,n.obs=n.obs,communality=c(rep(0,n)),loadings=matrix(rep(0,n*n),ncol=n),fit=0,fit.off=0)} else { result <- list(values=c(rep(0,n)),rotation=rotate,n.obs=n.obs,communality=c(rep(0,n)),loadings=matrix(rep(0,n*n),ncol=n),residual=matrix(rep(0,n*n),ncol=n),fit=0,fit.off=0)} eigens <- eigen(r) #call the eigen value decomposition routine result$values <- eigens$values eigens$values[ eigens$values < .Machine$double.eps] <- .Machine$double.eps #added May 14, 2009 to fix case of singular matrices loadings <- eigens$vectors %*% sqrt(diag(eigens$values)) if(nfactors >0) {loadings <- loadings[,1:nfactors]} else {nfactors <- n} if (nfactors>1) {communalities <- rowSums(loadings^2)} else {communalities <- loadings^2 } names(communalities) <- colnames(r) # 2009.02.10 Make sure this is a named vector -- correction by Gumundur Arnkelsson #added January 19, 2009 to flip based upon colSums of loadings if (nfactors >1) {sign.tot <- vector(mode="numeric",length=nfactors) sign.tot <- sign(colSums(loadings)) sign.tot[sign.tot==0] <- 1 loadings <- loadings %*% diag(sign.tot) } else { if (sum(loadings) <0) {loadings <- -as.matrix(loadings)} else {loadings <- as.matrix(loadings)} colnames(loadings) <- "PC1" } colnames(loadings) <- paste("PC",1:nfactors,sep='') rownames(loadings) <- rownames(r) Phi <- NULL if(rotate != "none") {if (nfactors > 1) { if (rotate=="varimax" | rotate== "Varimax" |rotate=="quartimax") { rotated <- do.call(rotate,list(loadings)) loadings <- rotated$loadings colnames(loadings) <- paste("RC",1:nfactors,sep='') Phi <- NULL} else { if ((rotate=="promax")|(rotate=="Promax")) {pro <- Promax(loadings) loadings <- pro$loadings Phi <- pro$Phi} else { if (rotate == "cluster") {loadings <- Varimax(loadings)$loadings pro <- target.rot(loadings) loadings <- pro$loadings colnames(loadings) <- paste("TC",1:nfactors,sep='') Phi <- pro$Phi} else { if (rotate =="oblimin"| rotate=="quartimin" | rotate== "simplimax") { if (!require(GPArotation)) {warning("I am sorry, to do these rotations requires the GPArotation package to be installed") Phi <- NULL} else { ob <- do.call(rotate,list(loadings) ) loadings <- ob$loadings colnames(loadings) <- paste("TC",1:nfactors,sep='') Phi <- ob$Phi} } }}} }} #just in case the rotation changes the order of the components, sort them by size of eigen value if(nfactors >1) { ev.rotated <- diag(t(loadings) %*% loadings) ev.order <- order(ev.rotated,decreasing=TRUE) loadings <- loadings[,ev.order]} if(!is.null(Phi)) {Phi <- Phi[ev.order,ev.order] } #January 20, 2009 but, then, we also need to change the order of the rotation matrix! signed <- sign(colSums(loadings)) c.names <- colnames(loadings) signed[signed==0] <- 1 loadings <- loadings %*% diag(signed) #flips factors to be in positive direction but loses the colnames colnames(loadings) <- c.names if(!is.null(Phi)) {Phi <- diag(signed) %*% Phi %*% diag(signed) colnames(Phi) <- rownames(Phi) <- c.names} class(loadings) <- "loadings" result$n.obs <- n.obs stats <- factor.stats(r,loadings,Phi,n.obs) class(result) <- c("psych", "principal") result$fn <- "principal" result$loadings <- loadings result$Phi <- Phi result$Call <- cl result$communality <- communalities #result$stats <- stats result$R2 <- stats$R2 result$objective <- stats$objective result$residual <- stats$residual result$fit <- stats$fit result$fit.off <- stats$fit.off result$factors <- stats$factors result$dof <- stats$dof result$null.dof <- stats$null.dof result$null.model <- stats$null.model result$criteria <- stats$criteria result$STATISTIC <- stats$STATISTIC result$PVAL <- stats$PVAL result$weights <- stats$weights result$r.scores <- stats$r.scores if(scores) {result$scores <- factor.scores(scale(x.matrix),loadings) } return(result) } # modified August 10, 2007 # modified Feb 2, 2008 to calculate scores and become a factanal class #Modified June 8,2008 to get chi square values to work or default to statistic if n.obs==NULL. #modified Jan 2, 2009 to report the correlations between oblique factors #modified December 30 to let n.obs ==NA rather than NULL to be compatible with factanal #modified Jan 14, 2009 to change phi to Phi to avoid confusion #modified Jan 19, 2009 to allow for promax rotations #modified May 15, 2009 to allow for singular matrices #correct August 25, 2009 to return result$values #modified August 25 to return result$stats