<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-component-method-part-two/"></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-component-method-part-two/"> R – Aaron Schlegel</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>In the <a target="_blank" rel="nofollow" href="http://wp.me/p4aZEo-5Qv">first post on factor analysis</a>, we examined computing the estimated covariance matrix <span data-display="false" class="wp-katex-eq">S</span> of the rootstock data and proceeded to find two factors that fit most of the variance of the data using the principal component method. However, the variables in the data are not on the same scale of measurement, which can cause variables with comparatively large variances to dominate the diagonal of the covariance matrix and the resulting factors. The correlation matrix, therefore, makes more intuitive sense to employ in factor analysis. In fact, as we saw previously, most packages available in R default to using the correlation matrix when performing factor analysis. There are several benefits to using <span data-display="false" class="wp-katex-eq">R</span> over <span data-display="false" class="wp-katex-eq">S</span>, not only that it scales non-commensurate variables, but it is also easier to calculate the factors as the matrix does not need to be decomposed and estimated like <span data-display="false" class="wp-katex-eq">S</span>.</p>
<h6>Factor Analysis with the Correlation Matrix</h6>
<p>Similar to factor analysis with the covariance matrix, we estimate <span data-display="false" class="wp-katex-eq">\Lambda</span> which is <span data-display="false" class="wp-katex-eq">p \times m</span> where <span data-display="false" class="wp-katex-eq">D</span> is a diagonal matrix of the <span data-display="false" class="wp-katex-eq">m</span> largest eigenvalues of <span data-display="false" class="wp-katex-eq">R</span>, and <span data-display="false" class="wp-katex-eq">C</span> is a matrix of the corresponding eigenvectors as columns.</p>
<p><span data-display="true" class="wp-katex-eq"> \hat{\Lambda} = CD^{1/2} = (\sqrt{\theta_1}c_1, \sqrt{\theta_2}c_2, \cdots, \sqrt{\theta_m}c_m)</span></p>
<p>Where <span data-display="false" class="wp-katex-eq">\theta_1, \theta_2, \cdots, \theta_m</span> are the largest eigenvalues of <span data-display="false" class="wp-katex-eq">R</span>.</p>
<p>Thus the correlation matrix <span data-display="false" class="wp-katex-eq">R</span> does not require decomposition, and we can proceed directly to finding the eigenvalues and eigenvectors of <span data-display="false" class="wp-katex-eq">R</span>.</p>
<p>Load the rootstock data and name the columns. From the previous post:</p>
<p>The rootstock data contains growth measurements of six different apple tree rootstocks from 1918 to 1934 (Andrews and Herzberg 1985, pp. 357-360) and 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>
<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>Compute the correlation matrix of the data.</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>Then find the eigenvalues and eigenvectors of <span data-display="false" class="wp-katex-eq">R</span>.</p>
<pre>
r.eigen &lt;- eigen(R)
r.eigen
## $values
## [1] 2.78462702 1.05412174 0.11733950 0.04391174
## 
## $vectors
##           [,1]       [,2]       [,3]       [,4]
## [1,] 0.4713465  0.5600120  0.6431731  0.2248274
## [2,] 0.5089667  0.4544775 -0.7142114 -0.1559013
## [3,] 0.5243109 -0.4431448  0.2413716 -0.6859012
## [4,] 0.4938456 -0.5324091 -0.1340527  0.6743048
</pre>
<p>We can check the proportion of each eigenvalue respective to the total sum of the eigenvalues.</p>
<p><span data-display="true" class="wp-katex-eq"> \frac{\sum^p_{i=1} \hat{\lambda}^2_{ij}}{tr(R)} = \frac{\theta_j}{p}</span></p>
<p>Where <span data-display="false" class="wp-katex-eq">p</span> is the number of variables. The quick and dirty loop below finds the proportion of the total for each eigenvalue and the cumulative proportion.</p>
<pre>
cumulative.proportion &lt;- 0
prop &lt;- c()
cumulative &lt;- c()
for (i in r.eigen$values) {
  proportion &lt;- i / dim(root[,2:5])[2]
  cumulative.proportion &lt;- cumulative.proportion + proportion
  
  prop &lt;- append(prop, proportion)
  cumulative &lt;- append(cumulative, cumulative.proportion)
}
data.frame(cbind(prop, cumulative))
##         prop cumulative
## 1 0.69615676  0.6961568
## 2 0.26353043  0.9596872
## 3 0.02933488  0.9890221
## 4 0.01097793  1.0000000
</pre>
<p>As in the case of the covariance matrix, the first two factors account for nearly all of the sample variance and thus can proceed with <span data-display="false" class="wp-katex-eq">m = 2</span> factors.</p>
<p>The eigenvectors corresponding to the two largest eigenvalues are multiplied by the square roots of their respective eigenvalues as seen earlier to obtain the factor loadings.</p>
<pre>
factors &lt;- t(t(r.eigen$vectors[,1:2]) * sqrt(r.eigen$values[1:2]))
round(factors, 2)
##      [,1]  [,2]
## [1,] 0.79  0.57
## [2,] 0.85  0.47
## [3,] 0.87 -0.45
## [4,] 0.82 -0.55
</pre>
<p>Computing the communality remains the same as in the covariance setting.</p>
<pre>
h2 &lt;- rowSums(factors^2)
</pre>
<p>The specific variance when factoring <span data-display="false" class="wp-katex-eq">R</span> is <span data-display="false" class="wp-katex-eq">1 – \hat{h}^2_i</span>.</p>
<pre>
u2 &lt;- 1 - h2
</pre>
<p>According to the documentation of the <code>principal()</code> function (called by `?principal), there is another statistic called complexity, which is the number of factors on which a variable has moderate or high loadings (Rencher, 2002 pp. 431), that is found by:</p>
<p><span data-display="true" class="wp-katex-eq"> \frac{(\sum^m_{i=1} \hat{\lambda}^2_i)^2}{\sum^m_{i=1} \hat{\lambda}_i^4}</span></p>
<p>In the most simple structure, the complexity of all the variables is <span data-display="false" class="wp-katex-eq">1</span>. The complexity of the variables is reduced by performing rotation which will be seen later.</p>
<pre>
com &lt;- rowSums(factors^2)^2 / rowSums(factors^4)
com
## [1] 1.831343 1.553265 1.503984 1.737242
</pre>
<pre>
mean(com)
## [1] 1.656459
</pre>
<p>As seen in the previous post, the <code>principal()</code> function from the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/psych/">psych package</a> performs factor analysis with the principal component method.</p>
<pre>
library(psych)
</pre>
<p>Since we are using <span data-display="false" class="wp-katex-eq">R</span> instead of <span data-display="false" class="wp-katex-eq">S</span>, the <code>covar</code> argument remains <code>FALSE</code> by default. No rotation is done for now, so the <code>rotate</code> argument is set to <code>none</code>.</p>
<pre>
root.fa &lt;- principal(root[,2:5], nfactors = 2, rotate = 'none')
root.fa
## Principal Components Analysis
## Call: principal(r = root[, 2:5], nfactors = 2, rotate = &quot;none&quot;)
## Standardized loadings (pattern matrix) based upon correlation matrix
##                               PC1   PC2   h2    u2 com
## Trunk.Girth.4.Years          0.79  0.57 0.95 0.051 1.8
## Ext.Growth.4.Years           0.85  0.47 0.94 0.061 1.6
## Trunk.Girth.15.Years         0.87 -0.45 0.97 0.027 1.5
## Weight.Above.Ground.15.Years 0.82 -0.55 0.98 0.022 1.7
## 
##                        PC1  PC2
## SS loadings           2.78 1.05
## Proportion Var        0.70 0.26
## Cumulative Var        0.70 0.96
## Proportion Explained  0.73 0.27
## Cumulative Proportion 0.73 1.00
## 
## Mean item complexity =  1.7
## 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>The output of the <code>principal()</code> function agrees with our calculations.</p>
<h6>Factor Rotation with Varimax Rotation</h6>
<p>Rotation moves the axes of the loadings to produce a more simplified structure of the factors to improve interpretation. Therefore the goal of rotation is to find an interpretable pattern of the loadings where variables are clustered into groups corresponding to the factors. We will see that a successful rotation yields a complexity closer to <span data-display="false" class="wp-katex-eq">1</span>, which denotes the variables load highly on only one factor.</p>
<p>One of the most common approaches to rotation is <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Varimax_rotation">varimax rotation</a>, which is a type of orthogonal rotation (axes remain perpendicular). The varimax technique seeks loadings that maximize the variance of the squared loadings in each column of the rotated matrix <span data-display="false" class="wp-katex-eq">\hat{\Lambda}*</span>.</p>
<p>The <code>varimax()</code> function is used to find the rotated factor loadings. For those interested, the R code for the <code>varimax()</code> function can be found <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Talk:Varimax_rotation">here</a>.</p>
<pre>
factors.v &lt;- varimax(factors)$loadings
round(factors.v, 2)
## 
## Loadings:
##      [,1] [,2]
## [1,] 0.16 0.96
## [2,] 0.28 0.93
## [3,] 0.94 0.29
## [4,] 0.97 0.19
## 
##                 [,1]  [,2]
## SS loadings    1.928 1.907
## Proportion Var 0.482 0.477
## Cumulative Var 0.482 0.959
</pre>
<p>The varimax rotation was rather successful in finding a rotation that simplified the complexity of the variables. The first two variables now load highly on the second factor while the remaining two variables load primarily on the first factor.</p>
<p>Since we used an orthogonal rotation technique, the communalities will not change.</p>
<pre>
h2.v &lt;- rowSums(factors.v^2)
h2.v
## [1] 0.9492403 0.9390781 0.9725050 0.9779253
</pre>
<pre>
h2
## [1] 0.9492403 0.9390781 0.9725050 0.9779253
</pre>
<p>Thus the specific variances will also be unchanged.</p>
<pre>
u2.v &lt;- 1 - h2.v
u2.v
## [1] 0.05075965 0.06092192 0.02749496 0.02207470
</pre>
<pre>
u2
## [1] 0.05075965 0.06092192 0.02749496 0.02207470
</pre>
<p>As stated previously, the complexity of the variables on the rotated factors should be closer to <span data-display="false" class="wp-katex-eq">1</span> compared to the non-rotated complexity.</p>
<pre>
com.v &lt;- rowSums(factors.v^2)^2 / rowSums(factors.v^4)
com.v
## [1] 1.054355 1.179631 1.185165 1.074226
</pre>
<pre>
mean(com.v)
## [1] 1.123344
</pre>
<p>The complexity is rather close to <span data-display="false" class="wp-katex-eq">1</span> which provides us further acknowledgment the factors are now in a more simplified structure.</p>
<p>Setting the <code>rotation</code> argument to <code>varimax</code> in the <code>principal()</code> function outputs the rotated factors and corresponding statistics.</p>
<pre>
root.fa2 &lt;- principal(root[,2:5], nfactors = 2, rotation = 'varimax')
root.fa2
## Principal Components Analysis
## Call: principal(r = root[, 2:5], nfactors = 2, rotation = &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>
<h6>Interpretation of Factors</h6>
<p>The factor analysis performed on the rootstock data yielded two latent variables that fit and explain the variance of the data quite sufficiently. We see both variables relating to measurements at four years load heavily on factor 2 while the 15-year measurements load mainly on the first factor. Thus we could designate names for the factors, or latent variables, such as ‘15 years growth’ and ‘4 years growth’, respectively. There isn’t any standard way of ‘naming’ factors as the interpretation can vary widely between each case. In this example, the factors make intuitive sense based on how they load on the variables; however, factors resulting from a factor analysis may not always make logic sense to the original data. If the resulting factors do not seem logical, changes to the approach such as adjusting the number of factors or the threshold of the loadings deemed important, or even a different method of rotation can be done to improve interpretation.</p>
<h6>References</h6>
<p>Rencher, A. (2002). Methods of Multivariate Analysis (2nd ed.). Brigham Young University: John Wiley &amp; Sons, Inc.</p>
<p><a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Talk:Varimax_rotation">https://en.wikipedia.org/wiki/Talk:Varimax_rotation</a></p>
<p><a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Varimax_rotation">https://en.wikipedia.org/wiki/Varimax_rotation</a></p>
<p><a target="_blank" rel="nofollow" href="http://web.stanford.edu/class/psych253/tutorials/FactorAnalysis.html">http://web.stanford.edu/class/psych253/tutorials/FactorAnalysis.html</a></p>
<p>The post <a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/factor-analysis-principal-component-method-part-two/">Factor Analysis with the Principal Component Method Part Two</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-component-method-part-two/"> 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/BGtj0_NCuyU"></a></div>