exploratory factor analysis, oblique rotation, variance explained

Question

The question how to compute the variance explained by a factor model obtained through exploratory factors analysis pops up from time to time. A summary with many possibilities is here: Calculating variance explained by factors after exploratory factor analysis with oblique rotation in R

I am using R, with the psych package as shown by the code snippet:

library(psych)

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

The mean communalities amount: 0.7418824 The summary table also reports the 0.74 as the explained variance for the oblimin rotation:

Oblimin Rotation:

[...]
                       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

If I do the same analysis with an orthogonal rotation, e.g. varimax, then I receive the same values for the total variance explained.

Varimax Rotation:
[...]

                       MR1  MR3  MR2  MR5  MR4
SS loadings           2.44 1.45 1.06 0.91 0.81
Proportion Var        0.27 0.16 0.12 0.10 0.09
Cumulative Var        0.27 0.43 0.55 0.65 0.74
Proportion Explained  0.37 0.22 0.16 0.14 0.12
Cumulative Proportion 0.37 0.58 0.74 0.88 1.00

In my understanding, an orthogonal rotation redistributes the loadings among the factors but the total variance explained remains the same for the whole system. However, I am unsure about an oblique rotation as the factors overlap in variance they explain. In the above example I have remarkable correlation among some of the factors but I receive the same total variance explained compared to varimax. Is the total explained variance also in this case the same or I do something wrong? I would need a feedback. Thanks!

Update:

I also looked at the output of fa_results$loadings in both cases. It displays not the same summary values for varimax as print(fa_results) but the output deviates in the case of oblimin rotation to a greater extent. Which values are the correct ones?

Oblimin: Deviation in comparison to print(fa_results)

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

Varimax: Deviation in comparison to print(fa_results).

Varimax Rotation:
>fa_results$loadings
[...]
                 MR1   MR3   MR2   MR5   MR4
SS loadings    2.444 1.450 1.064 0.913 0.806
Proportion Var 0.272 0.161 0.118 0.101 0.090
Cumulative Var 0.272 0.433 0.551 0.652 0.742

Answer 1

Your core question seems to be this:

In my understanding, an orthogonal rotation redistributes the loadings among the factors but the total variance explained remains the same for the whole system.

This is explained in Watkins (2021, p.88):

The unrotated and rotated factor solutions will explain the same amount of total variance (computed with eigenvalues) but rotation spreads that variance across the factors to improve interpretation and parsimony.

So this seems to align with the results (and your interpretation) of your output.

Reference

Watkins, M. W. (2021). A step-by-step guide to exploratory factor analysis with R and RStudio (1st ed.). Routledge. https://doi.org/10.4324/9781003120001

Answer 2

The update part of the question I posted on stackoverflow that was migrated to stackexchange but a full in-depth answer was posted by Prof. William Revelle, maintainer and writer of the psych package.

The solution can be read here: Discrepancies in explained variance depending on the print method