<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/iterated-principal-factor-method-of-factor-analysis-with-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.aaronschlegel.com/iterated-principal-factor-method-factor-analysis-r/"> R – Aaron Schlegel</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>The iterated principal factor method is an extension of the <a target="_blank" rel="nofollow" href="http://wp.me/p4aZEo-5QF">principal factor method</a> that seeks improved estimates of the communality. As seen in the previous post on the principal factor method, initial estimates of <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> are found to obtain <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span> from which the factors are computed. In the iterated principal factor method, the initial estimates of the communality are used to find new communality estimates from the loadings in <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span> with the following:</p>
<p><span data-display="true" class="wp-katex-eq"> \hat{h}^2_i = \sum^m_{j=1} \hat{\lambda}^2_{ij} </span></p>
<p>The values of <span data-display="false" class="wp-katex-eq">\hat{h}^2_i</span> are then substituted into the diagonal of <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> and a new value of <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span> is found. This iteration continues until the communality estimates converge, though sometimes convergence does not occur. Once the estimates converge, the eigenvalues and eigenvectors are calculated from the iterated <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 to arrive at the factor loadings.</p>
<h6>The Iterated Principal Factor Method in R</h6>
<p>The iterated principal factor method is demonstrated on the rootstock data as in the previous posts on factor analysis for consistency and comparison of the various approaches. The rootstock data contain four variables representing measurements in different units taken at four and fifteen years growth of six different rootstocks. 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>We proceed with the correlation matrix <span data-display="false" class="wp-katex-eq">R</span> as the variables in the data are not measured commensurately. Using <span data-display="false" class="wp-katex-eq">R</span> over <span data-display="false" class="wp-katex-eq">S</span> is generally the preferred approach and is usually the default in most implementations (such as the <code>psych</code> package).</p>
<p>Calculate the correlation matrix.</p>
<pre>
R &lt;- cor(root[,2:5])
</pre>
<p>The initial estimates of the communality are found by computing the squared multiple correlation, which in the case of <span data-display="false" class="wp-katex-eq">R – \hat{\Psi}</span> is equal 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>.</p>
<pre>
R.smc &lt;- (1 - 1 / diag(solve(R)))
</pre>
<p>The estimates then replace the diagonal of <span data-display="false" class="wp-katex-eq">R</span>.</p>
<pre>
diag(R) &lt;- R.smc
</pre>
<p>The threshold for convergence of the communality is set to <span data-display="false" class="wp-katex-eq">0.001</span>. This error threshold is also the default in the <code>psych</code> package implementation of the iterated principal factor method. The <code>com.iter</code> object will be used to store the communality iterations.</p>
<pre>
min.error &lt;- .001
com.iter &lt;- c()
</pre>
<p><span data-display="false" class="wp-katex-eq">\hat{h}^2_i</span> is then found from <span data-display="false" class="wp-katex-eq">\hat{h}^2_i = \sum^m_{j=1} \hat{\lambda}^2_{ij}</span>, which is simply the sum of the diagonal of <span data-display="false" class="wp-katex-eq">R</span> (<span data-display="false" class="wp-katex-eq">R</span> will be replaced with <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span> in the iteration).</p>
<pre>
h2 &lt;- sum(diag(R))
error &lt;- h2
</pre>
<p>The following loop implements the iterated principal factor method using <span data-display="false" class="wp-katex-eq">R</span> with the estimated communalities found earlier. While our communality estimate remains above the error threshold of <span data-display="false" class="wp-katex-eq">0.001</span>, the loop will continue to calculate new values of <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span> by replacing the previous estimates of communality with new ones.</p>
<pre>
while (error > min.error) {
  r.eigen &lt;- eigen(R) # Get the eigenvalues and eigenvectors of R
  
  # The lambda object is updated upon each iteration using new estimates of the communality
  lambda &lt;- as.matrix(r.eigen$vectors[,1:2]) %*% diag(sqrt(r.eigen$values[1:2]))
  
  # R - Psi is then found by multiplying the lambda matrix by its transpose
  r.mod &lt;- lambda %*% t(lambda)
  r.mod.diag &lt;- diag(r.mod) # The diagonal of R - Psi is the new communality estimate
  
  # The sum of the new estimate is taken and compared with the previous estimate. If the
  # difference is less than the error threshold the loop stops
  h2.new &lt;- sum(r.mod.diag) 
  error &lt;- abs(h2 - h2.new)
  
  # If the difference between the previous and new estimate is not below the threshold, replace
  # the new estimate with the previous
  h2 &lt;- h2.new
  
  # Store the iteration value (the sum of the estimate) and replace the diagonal of R with the
  # diagonal of R - Psi found previously
  com.iter &lt;- append(com.iter, h2.new)
  diag(R) &lt;- r.mod.diag
}
</pre>
<p>We now have the final <span data-display="false" class="wp-katex-eq">\hat{\Lambda}</span>. Find the communality, specific variances and complexity and collect them into a <code>data.frame</code></p>
<pre>
h2 &lt;- rowSums(lambda^2)
u2 &lt;- 1 - h2
com &lt;- rowSums(lambda^2)^2 / rowSums(lambda^4)

iter.fa.loadings &lt;- data.frame(cbind(round(lambda,2), round(h2, 2), round(u2, 3), round(com, 2)))
colnames(iter.fa.loadings) &lt;- c('Factor 1', 'Factor 2', 'h2', 'u2', 'com')
</pre>
<p>The proportion of variance explained by the factors is found by:</p>
<p><span data-display="true" class="wp-katex-eq"> \frac{\theta_j}{tr(R – \hat{\Psi})} = \frac{\theta_j}{\sum^p_{i=1} \theta_i} </span></p>
<p>Where <span data-display="false" class="wp-katex-eq">\theta_j</span> is the <span data-display="false" class="wp-katex-eq">j</span>th eigenvalue of <span data-display="false" class="wp-katex-eq">R – \hat{\Psi}</span>. The cumulative variance of the factors when factoring <span data-display="false" class="wp-katex-eq">R</span> is found by:</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. Calculate these values and store in a <code>data.frame</code>.</p>
<pre>
prop.var &lt;- r.eigen$values[1:2] / sum(diag(R))
var.cumulative &lt;- r.eigen$values / 4

factor.var &lt;- data.frame(rbind(round(prop.var[1:2], 2), round(var.cumulative[1:2], 2)))
rownames(factor.var) &lt;- c('Proportion Explained', 'Cumulative Variance')
colnames(factor.var) &lt;- c('Factor 1', 'Factor 2')
factor.var
##                      Factor 1 Factor 2
## Proportion Explained     0.74     0.26
## Cumulative Variance      0.68     0.24
</pre>
<p>The iterated principal factor method of factor analysis is complete, and we can now print the results!</p>
<pre>
iter.fa.res &lt;- list(iter.fa.loadings, factor.var)
iter.fa.res
## [[1]]
##   Factor 1 Factor 2   h2    u2  com
## 1    -0.76     0.56 0.89 0.109 1.84
## 2    -0.82     0.46 0.88 0.116 1.57
## 3    -0.88    -0.42 0.94 0.055 1.44
## 4    -0.83    -0.52 0.96 0.042 1.68
## 
## [[2]]
##                      Factor 1 Factor 2
## Proportion Explained     0.74     0.26
## Cumulative Variance      0.68     0.24
</pre>
<p>Interpretation of the factor loadings should be held until the factors are rotated. We will also compare the results of the iterated approach to the principal component method. Let’s compare our results with the output of the <code>psych</code> package to verify.</p>
<h6>Iterated Principal Factor Method with the <code>psych</code> Package</h6>
<p>The function <code>fa()</code> available in the <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/psych/">psych package</a> defaults to the iterated approach. We keep the <code>rotate</code> argument set to <code>none</code> for now and the <code>fm</code> argument to <code>pa</code> (principal axis, another term for principal factors).</p>
<pre>
library(psych)
</pre>
<pre>
root.cor.fa &lt;- fa(root[,2:5], nfactors = 2, rotate = 'none', fm = 'pa')
root.cor.fa
## Factor Analysis using method =  pa
## Call: fa(r = root[, 2:5], nfactors = 2, rotate = &quot;none&quot;, fm = &quot;pa&quot;)
## Standardized loadings (pattern matrix) based upon correlation matrix
##                               PA1   PA2   h2    u2 com
## Trunk.Girth.4.Years          0.76  0.56 0.89 0.109 1.8
## Ext.Growth.4.Years           0.82  0.46 0.88 0.116 1.6
## Trunk.Girth.15.Years         0.88 -0.42 0.94 0.055 1.4
## Weight.Above.Ground.15.Years 0.83 -0.52 0.96 0.042 1.7
## 
##                        PA1  PA2
## SS loadings           2.71 0.97
## Proportion Var        0.68 0.24
## Cumulative Var        0.68 0.92
## 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.06 
## 
## The root mean square of the residuals (RMSR) is  0.01 
## The df corrected root mean square of the residuals is  NA 
## 
## The harmonic number of observations is  48 with the empirical chi square  0.03  with prob &lt;  NA 
## The total number of observations was  48  with MLE Chi Square =  2.75  with prob &lt;  NA 
## 
## Tucker Lewis Index of factoring reliability =  1.128
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                 PA1  PA2
## Correlation of scores with factors             0.99 0.96
## Multiple R square of scores with factors       0.97 0.92
## Minimum correlation of possible factor scores  0.94 0.85
</pre>
<p>The output matches our results other than an arbitrary scaling of <span data-display="false" class="wp-katex-eq">-1</span> on the first factor in our calculations (notice this does not affect the communality or other computations as the loadings are squared). We can also see the <code>fa()</code> function had the same iterations as our loop using the <code>com.iter</code> object from earlier.</p>
<pre>
com.iter
## [1] 3.553558 3.615064 3.645859 3.661351 3.669220 3.673293 3.675480 3.676730
## [9] 3.677517
</pre>
<pre>
root.cor.fa$communality.iterations
## [1] 3.553558 3.615064 3.645859 3.661351 3.669220 3.673293 3.675480 3.676730
## [9] 3.677517
</pre>
<h6>Rotation of Factors</h6>
<p>The factors should be rotated so the variables load highly on one factor to better identify the groupings of the variables. Rotation also yields a simple structure of the data which is denoted by the complexity value calculated previously and improves interpretation of the factors.</p>
<p>The <code>varimax()</code> function can be used to rotate our computed factor loadings. Varimax rotation is a type of orthogonal rotation, in which the perpendicular axes remain perpendicular and the communality remains the same after rotation. Orthogonal rotations also result in uncorrelated factor loadings which can be useful for interpretation.</p>
<pre>
lambda.v &lt;- varimax(lambda)$loadings
lambda.v
## 
## Loadings:
##      [,1]   [,2]  
## [1,] -0.182  0.926
## [2,] -0.295  0.893
## [3,] -0.931  0.281
## [4,] -0.963  0.177
## 
##                 [,1]  [,2]
## SS loadings    1.912 1.765
## Proportion Var 0.478 0.441
## Cumulative Var 0.478 0.919
</pre>
<p>Otherwise, setting the <code>rotate</code> argument to <code>varimax</code> in the <code>fa()</code> function will perform varimax rotation.</p>
<h6>Comparison of Iterated Principal Factors and Principal Component Method</h6>
<p>We saw the non-iterated principal factor approach previously, and the <a target="_blank" rel="nofollow" href="http://wp.me/p4aZEo-5QF">principal component method</a> reported similar results; however, factor loadings from principal components loaded slightly higher on their respective variables and represented the more cumulative variance of the original data. Let’s see if the iterated method performs any better to the principal component method.</p>
<pre>
root.cor.fa &lt;- fa(root[,2:5], nfactors = 2, rotate = 'varimax', fm = 'pa')
root.cor.fa
## Factor Analysis using method =  pa
## Call: fa(r = root[, 2:5], nfactors = 2, rotate = &quot;varimax&quot;, fm = &quot;pa&quot;)
## Standardized loadings (pattern matrix) based upon correlation matrix
##                               PA1  PA2   h2    u2 com
## Trunk.Girth.4.Years          0.18 0.93 0.89 0.109 1.1
## Ext.Growth.4.Years           0.30 0.89 0.88 0.116 1.2
## Trunk.Girth.15.Years         0.93 0.28 0.94 0.055 1.2
## Weight.Above.Ground.15.Years 0.96 0.18 0.96 0.042 1.1
## 
##                        PA1  PA2
## SS loadings           1.91 1.77
## Proportion Var        0.48 0.44
## Cumulative Var        0.48 0.92
## Proportion Explained  0.52 0.48
## Cumulative Proportion 0.52 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.06 
## 
## The root mean square of the residuals (RMSR) is  0.01 
## The df corrected root mean square of the residuals is  NA 
## 
## The harmonic number of observations is  48 with the empirical chi square  0.03  with prob &lt;  NA 
## The total number of observations was  48  with MLE Chi Square =  2.75  with prob &lt;  NA 
## 
## Tucker Lewis Index of factoring reliability =  1.128
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                 PA1  PA2
## Correlation of scores with factors             0.98 0.96
## Multiple R square of scores with factors       0.97 0.93
## Minimum correlation of possible factor scores  0.93 0.86
</pre>
<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>Similar to the previous case with the non-iterated method, the principal component approach resulted in factors that loaded higher on their respective variables and represents slightly more cumulative variance of the data. The difference between the methods is rather small, yet one would be inclined to use the principal component method results.</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/iterated-principal-factor-method-factor-analysis-r/">Iterated Principal Factor Method of Factor Analysis with 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/iterated-principal-factor-method-factor-analysis-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/g_35ZQ3Qi-w"></a></div>