"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