Different results when performing PCA in R with princomp() and principal ()

1
votes

I tried to use princomp() and principal() to do PCA in R with data set USArressts. However, I got two different results for loadings/rotaion and scores.

First, I centered and normalised the original data frame so it is easier to compare the outputs.

library(psych)

trans_func <- function(x){
  x <- (x-mean(x))/sd(x)
  return(x)
}

A <- USArrests
USArrests <- apply(USArrests, 2, trans_func)

princompPCA <- princomp(USArrests, cor = TRUE)
principalPCA <- principal(USArrests, nfactors=4 , scores=TRUE, rotate = "none",scale=TRUE)

Then I got the results for the loadings and scores using the following commands:

princompPCA$loadings
principalPCA$loadings

Could you please help me to explain why there is a difference? and how can we interprete these results?

3
rpcapsych
In ?principle documentation, "Unlike princomp, this returns a subset of just the best nfactors." – chinsoon12
Further on, "The regression weights are found from the inverse of the correlation matrix times the component loadings. This has the result that the component scores are standard scores (mean=0, sd = 1) of the standardized input. A comparison to the scores from princomp shows this difference. princomp does not, by default, standardize the data matrix, nor are the components themselves standardized." – chinsoon12

3 Answers

3
votes

At the very end of the help document of ?principal:

"The eigen vectors are rescaled by the sqrt of the eigen values to produce the component loadings more typical in factor analysis."

So principal returns the scaled loadings. In fact, principal produces a factor model estimated by the principal component method.

1
votes

In 4 years, I would like to provide a more accurate answer to this question. I use iris data as an example.

data = iris[, 1:4]

First, do PCA by the eigen-decomposition

eigen_res = eigen(cov(data))
l = eigen_res$values
q = eigen_res$vectors

Then the eigenvector corresponding to the largest eigenvalue is the factor loadings

q[,1]

We can treat this as a reference or the correct answer. Now we check the results by different r functions. First, by function 'princomp'

res1 = princomp(data)
res1$loadings[,1]
# compare with 
q[,1]

No problem, this function actually just return the same results as 'eigen'. Now move to 'principal'

library(psych)
res2 = principal(data, nfactors=4, rotate="none")
# the loadings of the first PC is
res2$loadings[,1]
# compare it with the results by eigendecomposition
sqrt(l[1])*q[,1] # re-scale the eigen vector by sqrt of eigen value

You may find they are still different. The problem is the 'principal' function does eigendecomposition on the correlation matrix by default. Note: PCA is not invariant with rescaling the variables. If you modify the code as

res2 = principal(data, nfactors=4, rotate="none", cor="cov")
# the loadings of the first PC is
res2$loadings[,1]
# compare it with the results by eigendecomposition
sqrt(l[1])*q[,1] # re-scale the eigen vector by sqrt of eigen value

Now, you will get the same results as 'eigen' and 'princomp'.

Summarize:

  1. If you want to do PCA, you'd better apply 'princomp' function.
  2. PCA is a special case of the Factor model or a simplified version of the factor model. It is just equivalent to eigendecomposition.
  3. We can apply PCA to get an approximation of a factor model. It doesn't care about the specific factors, i.e. epsilons in a factor model. So, if you change the number of factors in your model, you will get the same estimations of the loadings. It is different from the maximum likelihood estimation.
  4. If you are estimating a factor model, you'd better use 'principal' function, since it provides more functions, like rotation, calculating the scores by different methods, and so on.
  5. Rescale the loadings of a PCA model doesn't affect the results too much. Since you still project the data onto the same optimal direction, i.e. maximize the variation in the resulting PC.
0
votes
ev <- eigen(R) # R is a correlation matrix of DATA
ev$vectors %*% diag(ev$values) %*% t(ev$vectors)

pc <- princomp(scale(DATA, center = F, scale = T),cor=TRUE) 
p <-principal(DATA, rotate="none")  

#eigen values
ev$values^0.5
pc$sdev
p$values^0.5

#eigen vectors - loadings
ev$vectors
pc$loadings
p$weights %*% diag(p$values^0.5)

pc$loading %*% diag(pc$sdev)
p$loadings 

#weights
ee <- diag(0,2)
for (j in 1:2) {
 for (i in 1:2) {
  ee[i,j] <- ev$vectors[i,j]/p$values[j]^0.5
 }
};ee 

#scores
s <- as.matrix(scale(DATA, center = T, scale = T)) %*% ev$vectors
scale(s)
p$scores
scale(pc$scores)