How to obtain unstandardized scores in factor analysis (FA)?

Question

When conducting factor analysis (FA), retained factors are standardised. I want to know if it is possible to obtain the factor scores in the same metric as the untandardized observed variables. This is the only answer I could find to my question, and a similar one here. However, in that example, there is only ONE factor. It seems in general that such "un-standardisation" is not normally performed, as there are no references to it at all.

My questions are:

1. Is it possible to unstandardize my factor scores when there is more than one factor?
2. What is the mathematical explanation for that (whether possible or not)?
3. Will it be dependent on the orthogonality of the factors? (i.e., what if I use an oblique rotation and factors are correlated?)

I'm attaching a similar code in R to the referenced question, only with factors = 3 and oblique rotation. Is this calculation valid?

v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6)
v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)
v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6)
v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)
v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5)
v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)
observed_matrix <- cbind(v1,v2,v3,v4,v5,v6)
fa = factanal(observed_matrix, factors = 3, scores = "regression",  rotation='promax')
FactorScores <- fa$scores
FactorLoadings <- fa$loadings
UnstandardizedFactorScores<-rowMeans(observed_matrix)+apply(observed_matrix,1,sd)*FactorScores 

Answer 1

An overview of methods to compute component and factor scores notices on the so called "standardized" factor scores,

The scores computed ... are scaled: they have variances equal to or close to 1 (standardized or near standardized) - not the true factor variances (which equal the sum of squared structure loadings).

and

In FA (not PCA), regressionally computed factor scores will appear not quite "standardized" - will have variances not 1, but equal to the $\frac {SS_{regr}}{(n-1)}$ of regressing these scores by the variables.

Please read that answer for more about what variance of factor score different methods produce. If the variance is not $1$, you have to z-standardize the scores "once more time" if you want their variance become exactly $1$.

To rescale a z-standardized (with st. dev. = 1) variable to have an arbitrary variance we simply multiply its values by the st. deviation we need, as you know.

So, if $F_z$ is z-standardized factor scores and $s^2$ is the true factor variance, then $F_z s$ yields the unstandardized factor scores you want. Factor variance $s^2$ equals the sum of squared loadings of that factor. So, to obtain unstandardized factor scores is straightforward and easy. The easiness is probably the reason why it is rarely described in literature. To cite the referred answer again:

So, when you need to supply factor scores with the true factor's variance, multiply the scores (having standardized them to st.dev. 1) by the sq. root of that variance.

When factors were modeled orthogonal, you don't have to worry about their correlations, and perform the above unstandardization of the scores independently for each factor, whatever their number.

Oblique factors. When factors are nonorthogonal, their variances (which are the sums of squared loadings in the structure matrix) partly overlap. So, the above method, being still valid, yields unstandardized scores which variances partly dublicate each other (the sum of variances is larger than the real multivariate variability the factors explain). A Footnote 3 here expresses:

The variance of a factor (or component) is the sum of its squared structure loadings S, since they are covariances/correlations between variables and (unit-scaled) factors. After oblique rotation, factors can get correlated, and so their variances intersect. Consequently, the sum of their variances, SS in S, exceeds the overall communality explained, SS in A. If you want to reckon after factor i only the unique "clean" portion of its variance, multiply its variance by $1-R_i^2$ of the factor's dependence on the other factors, the quantity known as anti-image. It is the reciprocal of the i-th diagonal element of $\bf C^{-1}$. The sum of the "clean" portions of the variances will be less than the overall communality explained.

Having obtained that "net-portion" estimate of the variance, unstandardize the factor scores to it. You will have each factor scores variate having only its unique, clean part of the variance, but the sum of variances will underestimate the total (multivariate) factor variability.

---

Addendum. You are asking, "to obtain the factor scores in the same metric as the unstandardized observed variables"

If your FA was performed on the covariances between the input variables, my above answer applies directly, because the factor variances (sums of squared loadings) are then on the same scale as the original unscaled (only centered) variables.

If your FA was performed on the correlations (which is a more common case), then my above answer pertains, of course, to the variances of the standardized input data, - each variable having variance $1$. So this is not their original&unequal variances. But the overall multivariate variability of the original and the standardized data relate as $k=\frac{original}{standardized} = \frac{trace(Cov)}{p}$, where trace(Cov) is the diagonal sum of the covariance matrix and p is the number of variables. Multiplying the above-obtained unstandardized factor scores by constant $\sqrt k$ rescales them to your original variability.

(Please note that FA performed on correlations and performed on covariances yield different results which cannot be "converted" one into another. One has to decide in advance if multivariate analysis such as clustering, FA or PCA, will be done on correlations or covariances.)