Discrepancies in explained variance depending on the print method

Question

While performing an exploratory factor analysis, I want to see the explained variance by the factors. It is straightforward and I printed the output of the fa() method. One can get the same summary table as the output of the $loadings, as shown below. The problem is that the explained variance deviates in the 2 tables. What is the reason for this, is it a bug in the psych package? Which output is the correct one?

library(psych)

fa_results <- fa(df_sq1sq9, nfactors = 5, rotate = "oblimin")
print(fa_results)
fa_results$loadings

Please consider that I am not talking about the factor loadings but the Cumulative var row in the below tables. The discrepancies appear in the case of orthogonal and also in oblique rotations, and the deviations of the two outputs are greater in the case of oblique rotations. Therefore, the example is shown with oblimin rotation.

Oblimin Rotation:
>print(fa_results)

[...]
                       MR1  MR3  MR5  MR2  MR4
SS loadings           2.45 1.50 1.07 1.05 0.60
Proportion Var        0.27 0.17 0.12 0.12 0.07
Cumulative Var        0.27 0.44 0.56 0.68 0.74
Proportion Explained  0.37 0.22 0.16 0.16 0.09
Cumulative Proportion 0.37 0.59 0.75 0.91 1.00

 With factor correlations of 
      MR1  MR3  MR5   MR2   MR4
MR1  1.00 0.59 0.64 -0.01  0.30
MR3  0.59 1.00 0.33  0.12  0.21
MR5  0.64 0.33 1.00  0.28  0.11
MR2 -0.01 0.12 0.28  1.00 -0.17
MR4  0.30 0.21 0.11 -0.17  1.00

Oblimin Rotation:
>fa_results$loadings
[...]
                 MR1   MR3   MR5   MR2   MR4
SS loadings    2.310 1.350 0.975 1.038 0.563
Proportion Var 0.257 0.150 0.108 0.115 0.063
Cumulative Var 0.257 0.407 0.515 0.630 0.693

The question I also raised as an update in https://stats.stackexchange.com/questions/657692/exploratory-factor-analysis-oblique-rotation-variance-explained

Answer 1

The problem is in the print function in core R for an object of class loading. It automatically finds the column sums of the squared loadings. This ignores the fact that the factors are correlated. Thus, this value will vary by rotation (which it should not, because the total amount of variance accounted for remains constant across rotations.

library(psych)
f5 <- fa(bfi[1:25],5)  # 5 factors default to oblimin
f5l <- f5$loadings    
tr(t(f5l) %*% f5l %*% f5$Phi)  #this adjusts for the correlations


tr(t(f5l)%*% f5l)   #this is the same sum of squares and is fine for orthogonal factors but give an incorrect answer for oblique factors.

Compare this to the unrotated solution

f5.none <- fa(bfi[1:25],5, rotate="none")
sum(f5.none$Vaccounted[1,])   #total variance accounted for matches the oblique above.

Then, although you did not ask, you might look at the communalities which in the orthogal solution are just the row sums of squares. In the oblique case, I find them by the diagonal of the reproduced correlation matrix:

h2 <- diag(f5l %*% f5$Phi %*% t(f5l))
sum(h2) 

#is the same as the total variance accounted for from above.