#pca is just an alias for the principal function
"pca" <- function(r,nfactors=1,residuals=FALSE,rotate="varimax",n.obs = NA, covar=FALSE,scores=TRUE,missing=FALSE,impute="median",oblique.scores=TRUE,method="regression",...) {
principal(r=r,nfactors=nfactors,residuals=residuals,rotate=rotate,n.obs=n.obs,covar=covar,scores=scores, missing=missing, impute=impute,oblique.scores=oblique.scores, method=method,...)
}
#to make principal components more similar in output to normal pca in other systems.
"principal" <-
function(r,nfactors=1,residuals=FALSE,rotate="varimax",n.obs = NA, covar=FALSE,scores=TRUE,missing=FALSE,impute="median",oblique.scores=TRUE,method="regression",...) {
cl <- match.call()
n <- dim(r)[2]
if (!isCorrelation(r)) {
raw <- TRUE
n.obs <- dim(r)[1]
if(scores) {x.matrix <- as.matrix(r) #matrices are required for the substitution to work
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]]}
}}
# 2011.12.21 added the covar option
if(!covar) {r <- cor(r,use="pairwise")} else r <- cov(r,use="pairwise") # if given a rectangular matrix, then find the correlations or covariances first
} else {
raw <- FALSE
if(!is.matrix(r)) { r <- as.matrix(r)}
sds <- sqrt(diag(r)) #convert covariance matrices to correlation matrices
if(!covar) 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)}
#added 24/4/15 to stop with bad data and give more meaningful help
if(any(is.na(r))) {
bad <- TRUE
tempr <-r
wcl <-NULL
while(bad) {
wc <- table(which(is.na(tempr), arr.ind=TRUE)) #find the correlations that are NA
wcl <- c(wcl,as.numeric(names(which(wc==max(wc)))))
tempr <- r[-wcl,-wcl]
if(any(is.na(tempr))) {bad <- TRUE} else {bad <- FALSE}
}
cat('\nLikely variables with missing values are ',colnames(r)[wcl],' \n')
stop("I am sorry: missing values (NAs) in the correlation matrix do not allow me to continue.\nPlease drop those variables and try again." )
}
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,nrow=length(eigens$values))) #added May 2, 2016 for the weird case of a single variable with covariance > 1
if(nfactors > 0) {loadings <- loadings[,1:nfactors]} else {nfactors <- n}
if (nfactors > 1) {communalities <- rowSums(loadings^2)} else {communalities <- loadings^2 }
uniquenesses <- diag(r) - communalities # 2011.12.21 uniqueness is now found if covar is true
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
rot.mat <- NULL
if(rotate != "none") {if (nfactors > 1) {
if (rotate=="varimax" |rotate=="Varimax" | rotate=="quartimax" | rotate =="bentlerT" | rotate =="geominT" | rotate =="targetT" | rotate =="bifactor" | rotate =="TargetT"|
rotate =="equamax"| rotate =="varimin"|rotate =="specialT" | rotate =="Promax" | rotate =="promax"| rotate =="cluster" |rotate == "biquartimin" |rotate == "TargetQ" |rotate =="specialQ" ) {
Phi <- NULL
colnames(loadings) <- paste("RC",1:nfactors,sep='') #for rotated component
switch(rotate, #The orthogonal cases for GPArotation + ones developed for psych
varimax = {rotated <- stats::varimax(loadings,...) #varimax is from stats, the others are from GPArotation
loadings <- rotated$loadings
rot.mat <- rotated$rotmat},
Varimax = {if (!requireNamespace('GPArotation')) {stop("I am sorry, to do this rotation requires the GPArotation package to be installed")}
#varimax is from the stats package, Varimax is from GPArotations
#rotated <- do.call(rotate,list(loadings,...))
#rotated <- do.call(getFromNamespace(rotate,'GPArotation'),list(loadings,...))
rotated <- GPArotation::Varimax(loadings,...)
loadings <- rotated$loadings
rot.mat <- t(solve(rotated$Th))} ,
quartimax = {if (!requireNamespace('GPArotation')) {stop("I am sorry, to do this rotation requires the GPArotation package to be installed")}
#rotated <- do.call(rotate,list(loadings))
rotated <- GPArotation::quartimax(loadings,...)
loadings <- rotated$loadings
rot.mat <- t(solve(rotated$Th))} ,
bentlerT = {if (!requireNamespace('GPArotation')) {stop("I am sorry, to do this rotation requires the GPArotation package to be installed")}
#rotated <- do.call(rotate,list(loadings,...))
rotated <- GPArotation::bentlerT(loadings,...)
loadings <- rotated$loadings
rot.mat <- t(solve(rotated$Th))} ,
geominT = {if (!requireNamespace('GPArotation')) {stop("I am sorry, to do this rotation requires the GPArotation package to be installed")}
#rotated <- do.call(rotate,list(loadings,...))
rotated <- GPArotation::geominT(loadings,...)
loadings <- rotated$loadings
rot.mat <- t(solve(rotated$Th))} ,
targetT = {if (!requireNamespace('GPArotation')) {stop("I am sorry, to do this rotation requires the GPArotation package to be installed")}
rotated <- GPArotation::targetT(loadings,Tmat=diag(ncol(loadings)),...)
loadings <- rotated$loadings
rot.mat <- t(solve(rotated$Th))} ,
bifactor = {rot <- bifactor(loadings,...)
loadings <- rot$loadings
rot.mat <- t(solve(rot$Th))},
TargetT = {if (!requireNamespace('GPArotation')) {stop("I am sorry, to do this rotation requires the GPArotation package to be installed")}
rot <- GPArotation::targetT(loadings,Tmat=diag(ncol(loadings)),...)
loadings <- rot$loadings
rot.mat <- t(solve(rot$Th))},
equamax = {rot <- equamax(loadings,...)
loadings <- rot$loadings
rot.mat <- t(solve(rot$Th))},
varimin = {rot <- varimin(loadings,...)
loadings <- rot$loadings
rot.mat <- t(solve(rot$Th))},
specialT = {rot <- specialT(loadings,...)
loadings <- rot$loadings
rot.mat <- t(solve(rot$Th))},
Promax = {pro <- Promax(loadings,...)
loadings <- pro$loadings
Phi <- pro$Phi
rot.mat <- pro$rotmat},
promax = {pro <- stats::promax(loadings,...) #from stats
loadings <- pro$loadings
rot.mat <- pro$rotmat
ui <- solve(rot.mat)
Phi <- cov2cor(ui %*% t(ui))},
cluster = {loadings <- varimax(loadings,...)$loadings
pro <- target.rot(loadings)
loadings <- pro$loadings
Phi <- pro$Phi
rot.mat <- pro$rotmat},
biquartimin = {ob <- biquartimin(loadings,...)
loadings <- ob$loadings
Phi <- ob$Phi
rot.mat <- t(solve(ob$Th))},
TargetQ = {ob <- TargetQ(loadings,...)
loadings <- ob$loadings
Phi <- ob$Phi
rot.mat <- t(solve(ob$Th))},
specialQ = {ob <- specialQ(loadings,...)
loadings <- ob$loadings
Phi <- ob$Phi
rot.mat <- t(solve(pro$Th))})
} else {
colnames(loadings) <- paste("TC",1:nfactors,sep='') #for transformed components
#The following oblique cases all use GPArotation
if (rotate =="oblimin"| rotate=="quartimin" | rotate== "simplimax" | rotate =="geominQ" | rotate =="bentlerQ" |rotate == "targetQ" ) {
if (!requireNamespace('GPArotation')) {warning("I am sorry, to do these rotations requires the GPArotation package to be installed")
Phi <- NULL} else {
ob <- try(do.call(getFromNamespace(rotate,'GPArotation'),list(loadings,...)))
if(class(ob)== as.character("try-error")) {warning("The requested transformaton failed, Promax was used instead as an oblique transformation")
ob <- Promax(loadings)}
loadings <- ob$loadings
Phi <- ob$Phi
rot.mat <- t(solve(ob$Th))}
} else {message("Specified rotation not found, rotate='none' used")
colnames(loadings) <- paste("PC",1:nfactors,sep='') } #not rotated
}
}
}
#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,fm="pc")
class(result) <- c("psych", "principal")
result$fn <- "principal"
result$loadings <- loadings
result$Phi <- Phi
result$Call <- cl
result$communality <- communalities
result$uniquenesses <- uniquenesses
result$complexity <- stats$complexity
#result$stats <- stats
result$chi <- stats$chi
result$EPVAL <- stats$EPVAL
result$R2 <- stats$R2
result$objective <- stats$objective
result$residual <- stats$residual
result$rms <- stats$rms
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
result$rot.mat <- rot.mat
if(!is.null(Phi) && oblique.scores) {
result$Structure <- loadings %*% Phi} else {result$Structure <- loadings
}
if(scores && raw) {
# if(covar) { if(!is.null(Phi)) {Phi <- Phi %*% diag(sqrt(diag(t(loadings) %*% loadings))) } }
result$weights <- solve(r,result$Structure)
result$scores <- scale(x.matrix,scale=!covar) %*% result$weights}
# result$scores <- factor.scores(scale(x.matrix,scale=!covar),result$Structure,Phi=Phi,method=method,rho=r) # method = method added Nov 20, 2011
# result$weights<- result$scores$weights
# result$scores <- result$scores$scores}
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
#modified November 20, 2011 to allow multiple scoring approaches to give component scores
#modified December 21, 2011 to allow for finding the pc of covariances
#modified June 19, 2012 to fix a missing value problem reported by Neil Stewart
#modified February 9, 2013 to label the type of rotation/transformation (as is documented but not implemented)
#modified February 25, 2013 to make scores=TRUE the default
#modified 1/16/14 to output the structure matrix
#modified 8/15 to add rot.mat to all rotations
#modified 9/3/15 to label rotated PCs as RC and transformed as TC -- perhaps this was dropped out while fixing rotations?
#added the pca function as an alias to principal 6/8/16
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.