Compute component scores from principal$loadings directly in R

Question

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?

Answer 1

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
quadrature points: 21 

Optimization:
Convergence: 0 
max(|grad|): 2.9e-05 
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.