<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/factor-analysis-with-the-principal-factor-method-and-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.aaronschlegel.com/factor-analysis-principal-factor-method-r/"> R – Aaron Schlegel</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>As discussed in a previous post on the <a target="_blank" rel="nofollow" href="http://wp.me/p4aZEo-5Qv">principal component method of factor analysis</a>, the <span data-display="false" class="wp-katex-eq">\hat{\Psi}</span> term in the estimated covariance matrix <span data-display="false" class="wp-katex-eq">S</span>, <span data-display="false" class="wp-katex-eq">S = \hat{\Lambda} \hat{\Lambda}’ + \hat{\Psi}</span>, was excluded and we proceeded directly to factoring <span data-display="false" class="wp-katex-eq">S</span> and <span data-display="false" class="wp-katex-eq">R</span>. The principal factor method of factor analysis (also called the principal axis method) finds an initial estimate of <span data-display="false" class="wp-katex-eq">\hat{\Psi}</span> and factors <span data-display="false" class="wp-katex-eq">S – \hat{\Psi}</span>, or <span data-display="false" class="wp-katex-eq">R – \hat{\Psi}</span> for the correlation matrix. Rearranging the estimated covariance and correlation matrices with the estimated <span data-display="false" class="wp-katex-eq">p \times m</span> <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span> matrix yields:</p>
<p><span data-display="true" class="wp-katex-eq"> S – \hat{\Psi} = \hat{\Lambda} \hat{\Lambda}^\prime </span><br>
<span data-display="true" class="wp-katex-eq"> R – \hat{\Psi} = \hat{\Lambda} \hat{\Lambda}^\prime </span></p>
<p>Therefore the principal factor method begins with eigenvalues and eigenvectors of <span data-display="false" class="wp-katex-eq">S – \hat{\Psi}</span> or <span data-display="false" class="wp-katex-eq">R – \hat{\Psi}</span>. <span data-display="false" class="wp-katex-eq">\hat{\Psi}</span> is a diagonal matrix of the <span data-display="false" class="wp-katex-eq">i</span>th communality. As in the principal component method, the <span data-display="false" class="wp-katex-eq">i</span>th communality, <span data-display="false" class="wp-katex-eq">\hat{h}^2_i</span>, is equal to <span data-display="false" class="wp-katex-eq">s_{ii} – \hat{\psi}_i</span> for <span data-display="false" class="wp-katex-eq">S – \hat{\Psi}</span> and <span data-display="false" class="wp-katex-eq">1 – \hat{\psi}_i</span> for <span data-display="false" class="wp-katex-eq">R – \hat{\Psi}</span>. The diagonal of <span data-display="false" class="wp-katex-eq">S</span> or <span data-display="false" class="wp-katex-eq">R</span> is replaced by their respective communalities in <span data-display="false" class="wp-katex-eq">\hat{\psi}_i</span> which gives us the following forms:</p>
<p><span data-display="true" class="wp-katex-eq"> S – \hat{\Psi} = \begin{bmatrix} \hat{h}^2_1 &amp;&amp; s_{12} &amp;&amp; \cdots &amp;&amp; s_{1p} \\\ s_{21} &amp;&amp; \hat{h}^2_2 &amp;&amp; \cdots &amp;&amp; s_{2p} \\\ \vdots &amp;&amp; \vdots &amp;&amp; &amp;&amp; \vdots \\\ s_{p1} &amp;&amp; s_{p2} &amp;&amp; \cdots &amp;&amp; \hat{h}^2_p \\\ \end{bmatrix}</span><br>
<span data-display="true" class="wp-katex-eq"> R – \hat{\Psi} = \begin{bmatrix} \hat{h}^2_1 &amp;&amp; r_{12} &amp;&amp; \cdots &amp;&amp; r_{1p} \\\ r_{21} &amp;&amp; \hat{h}^2_2 &amp;&amp; \cdots &amp;&amp; r_{2p} \\\ \vdots &amp;&amp; \vdots &amp;&amp; &amp;&amp; \vdots \\\ r_{p1} &amp;&amp; r_{p2} &amp;&amp; \cdots &amp;&amp; \hat{h}^2_p \\\ \end{bmatrix}</span></p>
<p>An initial estimate of the communalities is made using the squared multiple correlation between the observation vector <span data-display="false" class="wp-katex-eq">y_i</span> and the other <span data-display="false" class="wp-katex-eq">p – 1</span> variables. The squared multiple correlation in the case of <span data-display="false" class="wp-katex-eq">R – \hat{\Psi}</span> is equivalent to the following:</p>
<p><span data-display="true" class="wp-katex-eq"> \hat{h}^2_i = 1 – \frac{1}{r^{ii}} </span></p>
<p>Where <span data-display="false" class="wp-katex-eq">r^{ii}</span> is the <span data-display="false" class="wp-katex-eq">i</span>th diagonal element of <span data-display="false" class="wp-katex-eq">R^{-1}</span>. In the case of <span data-display="false" class="wp-katex-eq">S – \hat{\Psi}</span>, the above is multiplied by the variance of the respective variable.</p>
<p>The factor loadings are then calculated by finding the eigenvalues and eigenvectors of the <span data-display="false" class="wp-katex-eq">R – \hat{\Psi}</span> or <span data-display="false" class="wp-katex-eq">S – \hat{\Psi}</span> matrix.</p>
<h6>Factor Analysis with the Principal Factor Method in R</h6>
<p>We will perform factor analysis using the principal factor method on the rootstock data as done previously with the principal component method to see if the approaches differ significantly. The data were obtained from the <a target="_blank" rel="nofollow">companion FTP site</a> of the book Methods of Multivariate Analysis by Alvin Rencher. The data contains four dependent variables as follows:</p>
<ul>
<li>trunk girth at four years (mm <span data-display="false" class="wp-katex-eq">\times</span> 100)</li>
<li>extension growth at four years (m)</li>
<li>trunk girth at 15 years (mm <span data-display="false" class="wp-katex-eq">\times</span> 100)</li>
<li>weight of tree above ground at 15 years (lb <span data-display="false" class="wp-katex-eq">\times</span> 1000)</li>
</ul>
<p>Load the data and name the columns</p>
<pre>
root &lt;- read.table('ROOT.DAT', col.names = c('Tree.Number', 'Trunk.Girth.4.Years', 'Ext.Growth.4.Years', 'Trunk.Girth.15.Years', 'Weight.Above.Ground.15.Years'))
</pre>
<p>Since the variables of the rootstock data were measured on different scales, we will proceed with using the correlation matrix <span data-display="false" class="wp-katex-eq">R</span> to perform factor analysis.<br>
Find the correlation matrix.</p>
<pre>
R &lt;- cor(root[,2:5])
round(R, 2)
##                              Trunk.Girth.4.Years Ext.Growth.4.Years
## Trunk.Girth.4.Years                         1.00               0.88
## Ext.Growth.4.Years                          0.88               1.00
## Trunk.Girth.15.Years                        0.44               0.52
## Weight.Above.Ground.15.Years                0.33               0.45
##                              Trunk.Girth.15.Years
## Trunk.Girth.4.Years                          0.44
## Ext.Growth.4.Years                           0.52
## Trunk.Girth.15.Years                         1.00
## Weight.Above.Ground.15.Years                 0.95
##                              Weight.Above.Ground.15.Years
## Trunk.Girth.4.Years                                  0.33
## Ext.Growth.4.Years                                   0.45
## Trunk.Girth.15.Years                                 0.95
## Weight.Above.Ground.15.Years                         1.00
</pre>
<p>Calculate and replace the diagonal of <span data-display="false" class="wp-katex-eq">R</span> with the estimated communalities.</p>
<pre>
R.smc &lt;- (1 - 1 / diag(solve(R)))
diag(R) &lt;- R.smc
round(R, 2)
##                              Trunk.Girth.4.Years Ext.Growth.4.Years
## Trunk.Girth.4.Years                         0.80               0.88
## Ext.Growth.4.Years                          0.88               0.81
## Trunk.Girth.15.Years                        0.44               0.52
## Weight.Above.Ground.15.Years                0.33               0.45
##                              Trunk.Girth.15.Years
## Trunk.Girth.4.Years                          0.44
## Ext.Growth.4.Years                           0.52
## Trunk.Girth.15.Years                         0.91
## Weight.Above.Ground.15.Years                 0.95
##                              Weight.Above.Ground.15.Years
## Trunk.Girth.4.Years                                  0.33
## Ext.Growth.4.Years                                   0.45
## Trunk.Girth.15.Years                                 0.95
## Weight.Above.Ground.15.Years                         0.91
</pre>
<p>Now that we have an initial estimate of the communalities, we can find the eigenvalues and eigenvectors of the <span data-display="false" class="wp-katex-eq">R – \hat{\Psi}</span> matrix.</p>
<pre>
r.eigen &lt;- eigen(R)
r.eigen$values
## [1]  2.64598796  0.90756969 -0.03306996 -0.09008523
</pre>
<p>The <span data-display="false" class="wp-katex-eq">R</span> matrix is no longer positive semidefinite due to replacing the diagonal with the communalities. Thus there are a few small negative eigenvalues. Since negative eigenvalues cannot be used to estimate <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span> (due to taking the square root of the <span data-display="false" class="wp-katex-eq">D</span> matrix), we can proceed with <span data-display="false" class="wp-katex-eq">m = 2</span>.<br>
An interesting note is when negative eigenvalues exist, the cumulative proportion of variance calculated from the eigenvalues will exceed <span data-display="false" class="wp-katex-eq">1</span> and then decline back to <span data-display="false" class="wp-katex-eq">1</span> after considering the negative eigenvalues. The following loop demonstrates this phenomenon:</p>
<pre>
tot.prop &lt;- 0
for (i in r.eigen$values) {
  tot.prop &lt;- tot.prop + i / sum(r.eigen$values)
  print(tot.prop)
}
## [1] 0.7713346
## [1] 1.035901
## [1] 1.026261
## [1] 1
</pre>
<p>Obtain the factor loadings as before by multiplying the square root of the first two eigenvalues by their respective eigenvectors.</p>
<pre>
r.lambda &lt;- as.matrix(r.eigen$vectors[,1:2]) %*% diag(sqrt(r.eigen$values[1:2]))
r.lambda
##            [,1]       [,2]
## [1,] -0.7402973  0.5421762
## [2,] -0.8034091  0.4520610
## [3,] -0.8770380 -0.4050851
## [4,] -0.8266112 -0.4951379
</pre>
<p>The communalities, specific variances and complexity of the factor loadings can then be calculated.</p>
<pre>
r.h2 &lt;- rowSums(r.lambda^2)
r.u2 &lt;- 1 - r.h2
com &lt;- rowSums(r.lambda^2)^2 / rowSums(r.lambda^4)
</pre>
<p>Collect the results into a <code>data.frame</code>.</p>
<pre>
cor.pa &lt;- data.frame(cbind(round(r.lambda, 2), round(r.h2, 2), round(r.u2, 3), round(com, 1)))
colnames(cor.pa) &lt;- c('PA1', 'PA2', 'h2', 'u2', 'com')
cor.pa
##     PA1   PA2   h2    u2 com
## 1 -0.74  0.54 0.84 0.158 1.8
## 2 -0.80  0.45 0.85 0.150 1.6
## 3 -0.88 -0.41 0.93 0.067 1.4
## 4 -0.83 -0.50 0.93 0.072 1.6
</pre>
<h6>Principal Factor Method with the <code>psych</code> Package</h6>
<p>The <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/psych/">psych package</a> also performs factor analysis using the principal method with the <code>fa()</code> function.</p>
<pre>
library(psych)
</pre>
<p>The <code>fa()</code> function performs the iterated principal factor method by default, which, as the name implies, iterates the initial communality estimates with those of the resulting <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span> matrix until they converge. This approach to factor analysis will be demonstrated in a future post. Setting the <code>max.iter</code> argument to 1 will output the non-iterated principal factor method results.</p>
<pre>
root.cor.fa &lt;- fa(root[,2:5], nfactors = 2, rotate = 'none', fm = 'pa', max.iter = 1)
root.cor.fa
## Factor Analysis using method =  pa
## Call: fa(r = root[, 2:5], nfactors = 2, rotate = &quot;none&quot;, max.iter = 1, 
##     fm = &quot;pa&quot;)
## Standardized loadings (pattern matrix) based upon correlation matrix
##                               PA1   PA2   h2    u2 com
## Trunk.Girth.4.Years          0.74  0.54 0.84 0.158 1.8
## Ext.Growth.4.Years           0.80  0.45 0.85 0.150 1.6
## Trunk.Girth.15.Years         0.88 -0.41 0.93 0.067 1.4
## Weight.Above.Ground.15.Years 0.83 -0.50 0.93 0.072 1.6
## 
##                        PA1  PA2
## SS loadings           2.65 0.91
## Proportion Var        0.66 0.23
## Cumulative Var        0.66 0.89
## Proportion Explained  0.74 0.26
## Cumulative Proportion 0.74 1.00
## 
## Mean item complexity =  1.6
## Test of the hypothesis that 2 factors are sufficient.
## 
## The degrees of freedom for the null model are  6  and the objective function was  4.19 with Chi Square of  187.92
## The degrees of freedom for the model are -1  and the objective function was  0.17 
## 
## The root mean square of the residuals (RMSR) is  0.02 
## The df corrected root mean square of the residuals is  NA 
## 
## The harmonic number of observations is  48 with the empirical chi square  0.24  with prob &lt;  NA 
## The total number of observations was  48  with MLE Chi Square =  7.26  with prob &lt;  NA 
## 
## Tucker Lewis Index of factoring reliability =  1.281
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                 PA1  PA2
## Correlation of scores with factors             0.98 0.93
## Multiple R square of scores with factors       0.95 0.86
## Minimum correlation of possible factor scores  0.90 0.72
</pre>
<p>The results of the <code>fa()</code> function align to our own other than an arbitrary scaling of the first factor by <span data-display="false" class="wp-katex-eq">-1</span>.</p>
<h6>Rotation of Factor Loadings</h6>
<p>Rotate the factors using <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Varimax_rotation">varimax rotation</a> to improve interpretation. Varimax rotation of the loadings can be done with the <code>varimax()</code> function, or in the <code>fa()</code> function by setting the <code>rotate</code> argument to <code>varimax</code>.</p>
<pre>
root.cor.fa.v &lt;- fa(root[,2:5], nfactors = 2, rotate = 'varimax', fm = 'pa', max.iter = 1)
root.cor.fa.v
## Factor Analysis using method =  pa
## Call: fa(r = root[, 2:5], nfactors = 2, rotate = &quot;varimax&quot;, max.iter = 1, 
##     fm = &quot;pa&quot;)
## Standardized loadings (pattern matrix) based upon correlation matrix
##                               PA1  PA2   h2    u2 com
## Trunk.Girth.4.Years          0.19 0.90 0.84 0.158 1.1
## Ext.Growth.4.Years           0.30 0.87 0.85 0.150 1.2
## Trunk.Girth.15.Years         0.92 0.29 0.93 0.067 1.2
## Weight.Above.Ground.15.Years 0.95 0.18 0.93 0.072 1.1
## 
##                        PA1  PA2
## SS loadings           1.87 1.68
## Proportion Var        0.47 0.42
## Cumulative Var        0.47 0.89
## Proportion Explained  0.53 0.47
## Cumulative Proportion 0.53 1.00
## 
## Mean item complexity =  1.1
## Test of the hypothesis that 2 factors are sufficient.
## 
## The degrees of freedom for the null model are  6  and the objective function was  4.19 with Chi Square of  187.92
## The degrees of freedom for the model are -1  and the objective function was  0.17 
## 
## The root mean square of the residuals (RMSR) is  0.02 
## The df corrected root mean square of the residuals is  NA 
## 
## The harmonic number of observations is  48 with the empirical chi square  0.24  with prob &lt;  NA 
## The total number of observations was  48  with MLE Chi Square =  7.26  with prob &lt;  NA 
## 
## Tucker Lewis Index of factoring reliability =  1.281
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                 PA1  PA2
## Correlation of scores with factors             0.97 0.94
## Multiple R square of scores with factors       0.94 0.87
## Minimum correlation of possible factor scores  0.88 0.75
</pre>
<h6>Principal Component and Principal Factor Methods Comparison</h6>
<p>The principal factor method (and iterated principal factor method) will usually yield results close to the principal component method if either the correlations or the number of variables is large (Rencher, 2002, pp. 424).</p>
<p>Perform the principal component method of factor analysis and compare with the principal factor method.</p>
<pre>
root.cor.pa &lt;- principal(root[,2:5], nfactors = 2, rotate = 'varimax')
root.cor.pa
## Principal Components Analysis
## Call: principal(r = root[, 2:5], nfactors = 2, rotate = &quot;varimax&quot;)
## Standardized loadings (pattern matrix) based upon correlation matrix
##                               RC1  RC2   h2    u2 com
## Trunk.Girth.4.Years          0.16 0.96 0.95 0.051 1.1
## Ext.Growth.4.Years           0.28 0.93 0.94 0.061 1.2
## Trunk.Girth.15.Years         0.94 0.29 0.97 0.027 1.2
## Weight.Above.Ground.15.Years 0.97 0.19 0.98 0.022 1.1
## 
##                        RC1  RC2
## SS loadings           1.94 1.90
## Proportion Var        0.48 0.48
## Cumulative Var        0.48 0.96
## Proportion Explained  0.50 0.50
## Cumulative Proportion 0.50 1.00
## 
## Mean item complexity =  1.1
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.03 
##  with the empirical chi square  0.39  with prob &lt;  NA 
## 
## Fit based upon off diagonal values = 1
</pre>
<p>Both methods achieved a simple structure of the loadings following rotation. The loadings from each method are rather similar and don’t differ significantly, though the principal component method yielded factors that load more heavily on the variables which the factors hypothetically represent. However, the factors resulting from the principal component method explain 96% of the cumulative variance compared to 89% from the principal factor method. Though not a drastic difference, one is inclined to proceed with the principal component method in this case as the factors account for almost all of the variance in the variables.</p>
<h6>References</h6>
<p>Rencher, A. (2002). Methods of Multivariate Analysis (2nd ed.). Brigham Young University: John Wiley &amp; Sons, Inc.</p>
<p>The post <a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/factor-analysis-principal-factor-method-r/">Factor Analysis with the Principal Factor Method and R</a> appeared first on <a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/">Aaron Schlegel</a>.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.aaronschlegel.com/factor-analysis-principal-factor-method-r/"> R – Aaron Schlegel</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/hOYQlcZPG9Y"></a></div>