"principal" <- function(r,nfactors=1,residuals=FALSE,rotate="varimax",n.obs = NULL, scores=FALSE,missing=FALSE,impute="median", digits=2) { 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 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 <- round(eigens\$values,digits) loadings <- eigens\$vectors %*% sqrt(diag(eigens\$values)) if(nfactors >0) {loadings <- loadings[,1:nfactors]} else {nfactors <- n} if(nfactors >1) { maxabs <- apply(apply(loadings,2,abs),2,which.max) sign.max <- vector(mode="numeric",length=nfactors) for (i in 1: nfactors) {sign.max[i] <- sign(loadings[maxabs[i],i])} loadings <- loadings %*% diag(sign.max) #sign of largest loading is positive } else {if(nfactors >0) {mini <- min(loadings) maxi <- max(loadings) if (abs(mini) > maxi) {loadings <- -loadings } loadings <- as.matrix(loadings) } } colnames(loadings) <- paste("PC",1:nfactors,sep='') rownames(loadings) <- rownames(r) if(rotate != "none") {if (nfactors > 1) { if (rotate=="varimax") { loadings <- varimax(loadings)\$loadings } else { if (rotate=="promax") {loadings <- promax(loadings)\$loadings } else { if (rotate =="oblimin") { if (!require(GPArotation)) {warning("I am sorry, to do oblimin rotations requires the GPArotation package to be installed") phi <- NULL} else { ob <- oblimin(loadings) loadings <- ob\$loadings phi <- ob\$Phi} } }} }} class(loadings) <- "loadings" model <- loadings %*% t(loadings) result\$communality <- diag(model) result\$uniquenesses <- 1 - result\$communality residual<- r - model if (residuals) {result\$residual <-residual} diag(model) <- 1 model.inv <- solve(model) m.inv.r <- model.inv %*% r result\$dof <- n * (n-1)/2 - n * nfactors + (nfactors *(nfactors-1)/2) result\$objective <- sum(diag((m.inv.r))) - log(det(m.inv.r)) -n if (!is.null(n.obs)) {result\$STATISTIC <- result\$objective * (n.obs-1) -(2 * n + 5)/6 -(2*nfactors)/3 if (result\$dof >0) { result\$PVAL <- pchisq(result\$STATISTIC, result\$dof, lower.tail = FALSE)} else result\$PVAL < - NA} r2 <- sum(r*r) rstar2 <- sum(residual*residual) diag(residual) <- 0 rstar.off <- sum(residual*residual) r2.off <- r2 - n result\$fit <- round(1-rstar2/r2,digits) result\$fit.off <- round(1-rstar.off/r2.off,digits) result\$loadings <- round(loadings,digits) if (rotate =="oblimin") {result\$phi <- phi} result\$factors <- nfactors class(result) <- "factanal" if(scores) {result\$scores <- factor.scores(scale(x.matrix),loadings) } return(result) } # modified August 10, 2007 #last modified Feb 2, 2008 to calculate scores and become a factanal class