# Compute component scores from principal$loadings directly in R I am using a polychoric correlation matrix to run PCA, so I cannot obtain the scores directly from the function. I am currently manually plugging in the numbers from the principal$loadings matrix to calculate the factor scores, like this:

            RC1    RC2
Q1  0.323  0.808
Q2 -0.306  0.867
Q3 ......
.......

RC1 <- Q1*0.323-Q2*0.306+...
RC2 <- Q1*0.808+Q2*0.867+...


I think that there should be a better way to perform this calculation. Is there any way that will allow me to compute this easier?

• Providing a minimal working example (MWE) would be helpful. – Aleksandr Blekh Apr 4 '15 at 22:15
• Can't you just do matrix multiplication? Or am I missing something? – Aaron left Stack Overflow Apr 5 '15 at 0:55
• @Aaron, Arrgh I forgot about matrix multiplication!!! Thanks... – ceoec Apr 5 '15 at 5:19

If I correctly understood what you're doing, I think that the easier way of computing component scores is either one of the following approaches. The first and, likely, the easiest, one is to use item response theory (IRT)-based R package ltm (see this JSS paper for details and examples). The second approach is more traditional and uses popular psych package, along with hetcor package for computing polychoric correlations. Minimal reproducible examples (MRE) for both approaches can be found below.

Additionally, you can take a look at this answer and this answer, but keep in mind that fa.poly() function is deprecated, which IMHO seems to imply that hetcor() is now a preferred solution.

Approach 1:

library(ltm)
library(psych) # for 'bock' data set: http://www.inside-r.org/packages/cran/psych/docs/bock

data(bock)

lsat.fit <- rasch(lsat6)
summary(lsat.fit)
ltm::factor.scores(lsat.fit)


Results:

Call:
rasch(data = lsat6)

Model Summary:
log.Lik      AIC      BIC
-2466.938 4945.875 4975.322

Coefficients:
value std.err   z.vals
Dffclt.Q1 -3.6153  0.3266 -11.0680
Dffclt.Q2 -1.3224  0.1422  -9.3009
Dffclt.Q3 -0.3176  0.0977  -3.2518
Dffclt.Q4 -1.7301  0.1691 -10.2290
Dffclt.Q5 -2.7802  0.2510 -11.0743
Dscrmn     0.7551  0.0694  10.8757

Integration:
method: Gauss-Hermite

Optimization:
Convergence: 0
quasi-Newton: BFGS

Call:
rasch(data = lsat6)

Scoring Method: Empirical Bayes

Factor-Scores for observed response patterns:
Q1 Q2 Q3 Q4 Q5 Obs     Exp     z1 se.z1
1   0  0  0  0  0   3   2.364 -1.910 0.790
2   0  0  0  0  1   6   5.468 -1.439 0.793
3   0  0  0  1  0   2   2.474 -1.439 0.793
4   0  0  0  1  1  11   8.249 -0.959 0.801
5   0  0  1  0  0   1   0.852 -1.439 0.793
6   0  0  1  0  1   1   2.839 -0.959 0.801
7   0  0  1  1  0   3   1.285 -0.959 0.801
8   0  0  1  1  1   4   6.222 -0.466 0.816
9   0  1  0  0  0   1   1.819 -1.439 0.793
10  0  1  0  0  1   8   6.063 -0.959 0.801
11  0  1  0  1  1  16  13.288 -0.466 0.816
12  0  1  1  0  1   3   4.574 -0.466 0.816
13  0  1  1  1  0   2   2.070 -0.466 0.816
14  0  1  1  1  1  15  14.749  0.049 0.836
15  1  0  0  0  0  10  10.273 -1.439 0.793
16  1  0  0  0  1  29  34.249 -0.959 0.801
17  1  0  0  1  0  14  15.498 -0.959 0.801
18  1  0  0  1  1  81  75.060 -0.466 0.816
19  1  0  1  0  0   3   5.334 -0.959 0.801
20  1  0  1  0  1  28  25.834 -0.466 0.816
21  1  0  1  1  0  15  11.690 -0.466 0.816
22  1  0  1  1  1  80  83.310  0.049 0.836
23  1  1  0  0  0  16  11.391 -0.959 0.801
24  1  1  0  0  1  56  55.171 -0.466 0.816
25  1  1  0  1  0  21  24.965 -0.466 0.816
26  1  1  0  1  1 173 177.918  0.049 0.836
27  1  1  1  0  0  11   8.592 -0.466 0.816
28  1  1  1  0  1  61  61.235  0.049 0.836
29  1  1  1  1  0  28  27.709  0.049 0.836
30  1  1  1  1  1 298 295.767  0.593 0.862


Approach 2:

library(polycor)
library(psych)

data(bock)

# calculate polychoric correlations
#lsat.corr <- polychoric(lsat6, smooth=TRUE, global=TRUE, polycor=F, ML = FALSE, std.err=FALSE, progress=TRUE)
lsat.corr <- hetcor(lsat6, ML=TRUE)

# perform PCA
lsat.pca <- principal(r = lsat.corr$correlations, nfactors = 3, rotate = "Promax") summary(lsat.pca) # compute factor scores from the 'lsat6' data set with the 'lsat.pca' PCA solution lsat.pca$scores <- factor.scores(lsat6, lsat.pca)
print("")
print("Sample of calculated factor scores:")
print("")
print(head(lsat.pca$scores$scores))

# plot results
biplot.psych(lsat.pca)


Results:

Factor analysis with Call: principal(r = lsat.corr\$correlations, nfactors = 3, rotate = "Promax")

Test of the hypothesis that 3 factors are sufficient.
The degrees of freedom for the model is -2  and the objective function was  1.2

With component correlations of
PC1  PC3  PC2
PC1 1.00 0.09 0.00
PC3 0.09 1.00 0.15
PC2 0.00 0.15 1.00

Sample of calculated factor scores:

PC1        PC3       PC2
[1,] -2.656750 -2.6932246 -1.860959
[2,] -2.656750 -2.6932246 -1.860959
[3,] -2.656750 -2.6932246 -1.860959
[4,] -3.554983 -0.7170872 -1.031587
[5,] -3.554983 -0.7170872 -1.031587
[6,] -3.554983 -0.7170872 -1.031587


NOTE: In the second approach, I intentionally specified 3 factors to simplify the illustrative plot. However, I have also tried 5 factors and as of now still couldn't figure out why the results between the approaches are so different. Clarifications and/or explanations will be greatly appreciated.

• fa.poly is deprecated!? I have just learned it... Thank you for the details example! – ceoec Apr 5 '15 at 5:21
• I want to know why the results are so different too... – ceoec Apr 5 '15 at 5:27
• @ceoec: You're welcome! Thank you for upvoting/accepting my answer - I hope that it is helpful. In regard to fa.poly() being deprecated, that's the info I've found in the documentation, but it's just one reference, so I'm not 100% sure about that. In regard to results, I hope that other people will chime in and clarify that. – Aleksandr Blekh Apr 5 '15 at 5:36