3

This is my first post on this forum. I am new to R and factor analysis and I think I have a very simple question. I'm following Andy Field's "Discovering Statistics Using R" book (2012) for general guidance. Field recommends providing both the pattern and structure matrices when performing factor analysis using oblique rotation.

While the pattern matrix is simply the table of loadings, I am having more difficulty obtaining the structure matrix in R using the factanal() function.

To obtain the structure matrix for a PCA using the principal() function, Field provides the following formula: fit$loadings %*% fit$Phi, which multiplies the factor loading matrix by the factor correlation matrix. Although factanal() does store the factor correlations somewhere—as it provides them in a section of the general output (under loadings)—I cannot find an object called "Phi", "Factor Correlations", or an equivalent, within the output of the fit model (within R Studio). Thus, I don't know what term to put in Field's formula to replace Phi in order to get the structure matrix.

I have seen here and here that earlier, factanal() did not provide the factor correlations—yet now it does provide it in the output; I just don't know the right terminology to access it. Thanks for any help on this!

Edit: As per the book, I use the following formula with four factors and oblique rotation for PCA:

pc4 <- principal(raqData, nfactors = 4, rotate = "oblimin").

In that case when I double click on object "pc4" in R Studio there is an object called "Phi" and I have no trouble obtaining the structure matrix with pc4$loadings %*% pc4$Phi.

Next, I attempt to use EFA, instead of PCA, again with four factors and oblique rotation (promax). This step works and I can get the factor correlations with (cut some output for conciseness):

> fit <- factanal(mydata, 4, rotation="promax")
> fit
Call...
Uniquenesses...
Loadings...
SS loadings...
Factor Correlations:
        Factor1 Factor2 Factor3 Factor4
Factor1   1.000  -0.874   0.632   0.571
Factor2  -0.874   1.000  -0.118  -0.438
Factor3   0.632  -0.118   1.000   0.356
Factor4   0.571  -0.438   0.356   1.000
Test of the hypothesis that 4 factors are sufficient...

Next, I attempt the structure matrix formula and get the following error:

> fit$loadings %*% fit$Phi
Error in fit$loadings %*% fit$Phi : 
  requires numeric/complex matrix/vector arguments

But when I inspect the "fit" object obtained with factanal(), there is no object called 'Phi" (as there were for principal()). Not sure how to interpret the error above either. I have seen this error discussed here and here but I'm not clear about how to resolve it here.

rempsyc
  • 785
  • 5
  • 24
  • You might do better on an R site. –  Feb 16 '18 at 12:45
  • Thanks all for the info. I've seen on the help page cross-posting is not encouraged but the post can be migrated to another site, for example if it is better suited to the r section of stackoverflow per se. I suppose a moderator would need to do that for me? –  Feb 16 '18 at 14:34
  • There should be some text under your question which says "flag". If you do that and explain it then it could be migrated. I do not use StackOverflow myself so I cannot say if they would think it on-topic there but others will know. –  Feb 16 '18 at 14:51
  • Can you post code or a reproducible example? Did you use a rotation? If so, which one? – Jeremy Miles Feb 16 '18 at 17:30
  • Thanks so much for responding here Dr. Miles! I have added some more details in edit. Let me know if this is sufficient or if you need actual data. Best. – rempsyc Feb 16 '18 at 19:14
  • Hmmmm..... Try using the fa() function, in the psych package, instead. – Jeremy Miles Feb 16 '18 at 19:24
  • Your command would be: fit <- fa(mydata, 4, rotate="promax") – Jeremy Miles Feb 16 '18 at 19:25
  • Ok yep, this is what I ended up doing, but this confirms there's no way to call the correlation matrix from factanal()? Furthermore, as mean of cross-validation I attempted to compare the factor correlations with the two methods; I get value differences as big as .12... I know different methods often give different outcomes because they use different algorithms but this seems pretty big so I'm a bit confused! Changing the fm option of the fa() function to fm="ml" (as factanal uses) also changes the values but did not bring them closer to those given by factanal(). – rempsyc Feb 16 '18 at 19:47

1 Answers1

1

Though this is an old post but for helping other people out, the two components loadings and psi (not phi) in factanal output are fit$loadings and fit$uniqueness. You can access loadings matrix by writing fit$loadings[,1:q] where q is the factor you used. For example- I ran factanal on Harman23.cor data present in the R studio. And checked for loadings to get below output.

> Harman23.FA <- factanal(factors = 3, covmat = Harman23.cor)
> Harman23.FA$loadings

Loadings:
               Factor1 Factor2 Factor3
height          0.886   0.267  -0.130 
arm.span        0.937   0.195   0.280 
forearm         0.874   0.188         
lower.leg       0.877   0.230  -0.145 
weight          0.242   0.916  -0.106 
bitro.diameter  0.193   0.777         
chest.girth     0.137   0.755         
chest.width     0.261   0.646   0.159 

               Factor1 Factor2 Factor3
SS loadings      3.379   2.628   0.162
Proportion Var   0.422   0.329   0.020
Cumulative Var   0.422   0.751   0.771

Clearly $loadings is not returning a matrix. Instead you can use this:

> as.matrix(Harman23.FA$loadings[,1:3])
                 Factor1   Factor2      Factor3
height         0.8859687 0.2666293 -0.130079754
arm.span       0.9371891 0.1951310  0.280364451
forearm        0.8741374 0.1876181  0.089154498
lower.leg      0.8768750 0.2303982 -0.144815728
weight         0.2422827 0.9164637 -0.106482093
bitro.diameter 0.1925419 0.7766283 -0.020581101
chest.girth    0.1368022 0.7554074  0.003677235
chest.width    0.2605798 0.6458413  0.159110501

Note that the factor used here is 3 and hence while accessing loadings we have mentioned colums 1:3. Similarly for running factanal with a factor q, you should mention Harman23.FA$loadings[,1:q].

Also note that psi which is the uniquess matrix is a diagonal matrix. So you should access it like this

> psihat <- diag(Harman23.FA$uniquenesses)
> psihat
          [,1]  [,2]     [,3]      [,4]       [,5]      [,6]      [,7]      [,8]
[1,] 0.1270475 0.000 0.000000 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.005 0.000000 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
[3,] 0.0000000 0.000 0.192735 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
[4,] 0.0000000 0.000 0.000000 0.1570353 0.00000000 0.0000000 0.0000000 0.0000000
[5,] 0.0000000 0.000 0.000000 0.0000000 0.09005497 0.0000000 0.0000000 0.0000000
[6,] 0.0000000 0.000 0.000000 0.0000000 0.00000000 0.3593523 0.0000000 0.0000000
[7,] 0.0000000 0.000 0.000000 0.0000000 0.00000000 0.0000000 0.4106308 0.0000000
[8,] 0.0000000 0.000 0.000000 0.0000000 0.00000000 0.0000000 0.0000000 0.4896708

Otherwise your matrix multiplication fit$loadings %*% fit$Phi will be non-conformal

user11831949
  • 11
  • 3
  • 1
    As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers [in the help center](/help/how-to-ask). – Community Sep 16 '21 at 16:16