<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/principal-component-analysis-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://poissonisfish.wordpress.com/2017/01/23/principal-component-analysis-in-r/"> poissonisfish</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Principal component analysis (PCA) is routinely employed on a wide range of problems. From the detection of outliers to predictive modeling, PCA has the ability of projecting the observations described by <img title="p " alt="p " class="latex" src="https://s0.wp.com/latex.php?latex=p+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> variables into few orthogonal components defined at where the data ‘stretch’ the most, rendering a simplified overview. PCA is particularly powerful in dealing with multicollinearity and variables that outnumber the samples (<img title="p \gg n " alt="p \gg n " class="latex" src="https://s0.wp.com/latex.php?latex=p+%5Cgg+n+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">).</p>
<p><img data-image-description="" data-orig-file="https://i1.wp.com/poissonisfish.files.wordpress.com/2017/01/q7hip.gif?w=456&amp;ssl=1" data-orig-size="1000,400" data-image-title="q7hip" class="wp-image-835" data-comments-opened="1" data-recalc-dims="1" data-attachment-id="835" data-permalink="https://poissonisfish.wordpress.com/2017/01/23/principal-component-analysis-in-r/q7hip/" alt="q7hip" data-medium-file="https://i1.wp.com/poissonisfish.files.wordpress.com/2017/01/q7hip.gif?w=456&amp;ssl=1?w=300" data-large-file="https://i1.wp.com/poissonisfish.files.wordpress.com/2017/01/q7hip.gif?w=456&amp;ssl=1?w=1000" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i1.wp.com/poissonisfish.files.wordpress.com/2017/01/q7hip.gif?w=456&ssl=1"></p>
<p>It is an unsupervised method, meaning it will always look into the greatest sources of variation regardless of the data structure. Its counterpart, the partial least squares (PLS), is a supervised method and will perform the same sort of covariance decomposition, albeit building a user-defined number of components (frequently designated as latent variables) that minimize the SSE from predicting a specified outcome with an ordinary least squares (OLS). The PLS is worth an entire post and so I will refrain from casting a second spotlight.</p>
<p>In case PCA is entirely new to you, there is an excellent <a target="_blank" rel="nofollow" href="http://www.nature.com/nbt/journal/v26/n3/full/nbt0308-303.html">Primer from Nature Biotechnology</a> that I highly recommend. Notwithstanding the focus on life sciences, it should still be clear to others than biologists.</p>
<h3><strong>Mathematical foundation</strong></h3>
<p>There are numerous PCA formulations in the literature dating back as long as one century, but all in all PCA is pure linear algebra. One of the most popular methods is the singular value decomposition (SVD). The SVD algorithm breaks down a matrix <img title="X " alt="X " class="latex" src="https://s0.wp.com/latex.php?latex=X+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> of size <img title="n \times p " alt="n \times p " class="latex" src="https://s0.wp.com/latex.php?latex=n+%5Ctimes+p+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> into three pieces,</p>
<p><b><img title="X = U \Sigma V^T " alt="X = U \Sigma V^T " class="latex" src="https://s0.wp.com/latex.php?latex=X+%3D+U+%5CSigma+V%5ET+&bg=f9f9f9&%23038;fg=555555&%23038;s=2"></b></p>
<p>where <img title="U " alt="U " class="latex" src="https://s0.wp.com/latex.php?latex=U%C2%A0&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> is the matrix with the eigenvectors of <img title="X X^T " alt="X X^T " class="latex" src="https://s0.wp.com/latex.php?latex=X+X%5ET+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">, <img title="\Sigma " alt="\Sigma " class="latex" src="https://s0.wp.com/latex.php?latex=%5CSigma+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> is the diagonal matrix with the singular values and <img title="V^T " alt="V^T " class="latex" src="https://s0.wp.com/latex.php?latex=V%5ET%C2%A0&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> is the matrix with the eigenvectors of <img title="X^T X " alt="X^T X " class="latex" src="https://s0.wp.com/latex.php?latex=X%5ET+X+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">. These matrices are of size <img title="n \times n " alt="n \times n " class="latex" src="https://s0.wp.com/latex.php?latex=n+%5Ctimes+n+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">, <img title="n \times p " alt="n \times p " class="latex" src="https://s0.wp.com/latex.php?latex=n+%5Ctimes+p+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> and <img title="p \times p " alt="p \times p " class="latex" src="https://s0.wp.com/latex.php?latex=p+%5Ctimes+p+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">, respectively. The key difference of SVD compared to a matrix diagonalization (<img title="X = V \Sigma V^{-1} " alt="X = V \Sigma V^{-1} " class="latex" src="https://s0.wp.com/latex.php?latex=X+%3D+V+%5CSigma+V%5E%7B-1%7D+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">) is that <img title="U " alt="U " class="latex" src="https://s0.wp.com/latex.php?latex=U%C2%A0&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> and <img title="V^T " alt="V^T " class="latex" src="https://s0.wp.com/latex.php?latex=V%5ET%C2%A0&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> are distinct orthonormal (orthogonal and unit-vector) matrices.</p>
<p>PCA reduces the <img title="p " alt="p " class="latex" src="https://s0.wp.com/latex.php?latex=p+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> dimensions of your data set <img title="X " alt="X " class="latex" src="https://s0.wp.com/latex.php?latex=X+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> down to <img title="k " alt="k " class="latex" src="https://s0.wp.com/latex.php?latex=k+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> principal components (PCs). The scores from the first <img title="k " alt="k " class="latex" src="https://s0.wp.com/latex.php?latex=k+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> PCs result from multiplying the first <img title="k " alt="k " class="latex" src="https://s0.wp.com/latex.php?latex=k+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> columns of <img title="U " alt="U " class="latex" src="https://s0.wp.com/latex.php?latex=U+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> with the <img title="k \times k " alt="k \times k " class="latex" src="https://s0.wp.com/latex.php?latex=k+%5Ctimes+k+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> upper-left submatrix of <img title="\Sigma " alt="\Sigma " class="latex" src="https://s0.wp.com/latex.php?latex=%5CSigma+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">. The loading factors of the <img title="k^{th} " alt="k^{th} " class="latex" src="https://s0.wp.com/latex.php?latex=k%5E%7Bth%7D+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> PC are directly given in the <img title="k^{th} " alt="k^{th} " class="latex" src="https://s0.wp.com/latex.php?latex=k%5E%7Bth%7D+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> row in <img title="V^T " alt="V^T " class="latex" src="https://s0.wp.com/latex.php?latex=V%5ET+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">. Consequently, multiplying all scores and loadings recovers <img title="X " alt="X " class="latex" src="https://s0.wp.com/latex.php?latex=X+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">. You might as well keep in mind:</p>
<ul>
<li>PCs are ordered by the decreasing amount of variance explained</li>
<li>PCs are orthogonal <em>i.e.</em> uncorrelated to each other</li>
<li>The columns of <img title="X " alt="X " class="latex" src="https://s0.wp.com/latex.php?latex=X+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> should be mean-centered, so that the covariance matrix <img title="\approx X^{T} X " alt="\approx X^{T} X " class="latex" src="https://s0.wp.com/latex.php?latex=%5Capprox+X%5E%7BT%7D+X+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"></li>
<li>SVD-based PCA does not tolerate missing values (but there are solutions we will cover shortly)</li>
</ul>
<p>For a more elaborate explanation with introductory linear algebra, here is an excellent free <a target="_blank" rel="nofollow" href="https://datajobs.com/data-science-repo/SVD-Tutorial-%5BKirk-Baker%5D.pdf">SVD tutorial</a> I found online. At any rate, I guarantee you can master PCA without fully understanding the process.</p>
<h2><strong>Let’s get started with R</strong></h2>
<p>Although there is a plethora of PCA methods available for R, I will only introduce two,</p>
<ul>
<li><strong>prcomp</strong>, a default function from the R base package</li>
<li><strong>pcaMethods</strong>, a <a target="_blank" rel="nofollow" href="http://bioconductor.org/">Bioconductor</a> package that I frequently use for my own PCAs</li>
</ul>
<p>I will start by demonstrating that <strong>prcomp</strong> is based on the SVD algorithm, using the base <strong>svd</strong> function.</p>
<pre>
# Generate scaled 4*5 matrix with random std normal samples
set.seed(101)
mat &lt;- scale(matrix(rnorm(20), 4, 5))
dimnames(mat) &lt;- list(paste(&quot;Sample&quot;, 1:4), paste(&quot;Var&quot;, 1:5))

# Perform PCA
myPCA &lt;- prcomp(mat, scale. = F, center = F)
myPCA$rotation # loadings
myPCA$x # scores
</pre>
<p>By default, <strong>prcomp</strong> will retrieve <img title="min(n, p) " alt="min(n, p) " class="latex" src="https://s0.wp.com/latex.php?latex=min%28n%2C+p%29+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> PCs. Therefore, in our setting we expect having four PCs.The <strong>svd</strong> function will behave the same way:</p>
<pre>
# Perform SVD
mySVD &lt;- svd(mat)
mySVD # the diagonal of Sigma mySVD$d is given as a vector
sigma &lt;- matrix(0,4,4) # we have 4 PCs, no need for a 5th column
diag(sigma) &lt;- mySVD$d # sigma is now our true sigma matrix
</pre>
<p>Now that we have the PCA and SVD objects, let us compare the respective scores and loadings. We will compare the scores from the PCA with the product of <img title="U " alt="U " class="latex" src="https://s0.wp.com/latex.php?latex=U+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> and <img title="\Sigma " alt="\Sigma " class="latex" src="https://s0.wp.com/latex.php?latex=%5CSigma+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> from the SVD. In R, matrix multiplication is possible with the operator <strong>%*%</strong>. Next, we will directly compare the loadings from the PCA with <img title="V " alt="V " class="latex" src="https://s0.wp.com/latex.php?latex=V+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> from the SVD, and finally show that multiplying scores and loadings recovers <img title="X " alt="X " class="latex" src="https://s0.wp.com/latex.php?latex=X+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">. I rounded the results to five decimal digits since the results are not exactly the same! The function <strong>t </strong>retrieves a transposed matrix.</p>
<pre>
# Compare PCA scores with the SVD's U*Sigma
theoreticalScores &lt;- mySVD$u %*% sigma
all(round(myPCA$x,5) == round(theoreticalScores,5)) # TRUE

# Compare PCA loadings with the SVD's V
all(round(myPCA$rotation,5) == round(mySVD$v,5)) # TRUE

# Show that mat == U*Sigma*t(V)
recoverMatSVD &lt;- theoreticalScores %*% t(mySVD$v)
all(round(mat,5) == round(recoverMatSVD,5)) # TRUE

#Show that mat == scores*t(loadings)
recoverMatPCA &lt;- myPCA$x %*% t(myPCA$rotation)
all(round(mat,5) == round(recoverMatPCA,5)) # TRUE
</pre>
<h3></h3>
<h3>PCA of the wine data set</h3>
<p>Now that we established the association between SVD and PCA, we will perform PCA on real data. I found a wine data set at the <a target="_blank" rel="nofollow" href="http://archive.ics.uci.edu/ml/index.html">UCI Machine Learning Repository</a> that might serve as a good starting example.</p>
<pre>
wine &lt;- read.table(&quot;http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data&quot;, sep=&quot;,&quot;)
</pre>
<p>According to the <a target="_blank" rel="nofollow" href="http://archive.ics.uci.edu/ml/datasets/Wine">documentation</a>, these data consist of 13 physicochemical parameters measured in 178 wine samples from three distinct cultivars grown in Italy. Let’s check patterns in pairs of variables, and then see what a PCA does about that by plotting PC1 against PC2.</p>
<pre>
# Name the variables
colnames(wine) &lt;- c(&quot;Cvs&quot;,&quot;Alcohol&quot;,&quot;Malic acid&quot;,&quot;Ash&quot;,&quot;Alcalinity of ash&quot;, &quot;Magnesium&quot;, &quot;Total phenols&quot;, &quot;Flavanoids&quot;, &quot;Nonflavanoid phenols&quot;, &quot;Proanthocyanins&quot;, &quot;Color intensity&quot;, &quot;Hue&quot;, &quot;OD280/OD315 of diluted wines&quot;, &quot;Proline&quot;)

# The first column corresponds to the classes
wineClasses &lt;- factor(wine$Cvs)

# Use pairs
pairs(wine[,-1], col = wineClasses, upper.panel = NULL, pch = 16, cex = 0.5)
legend(&quot;topright&quot;, bty = &quot;n&quot;, legend = c(&quot;Cv1&quot;,&quot;Cv2&quot;,&quot;Cv3&quot;), pch = 16, col = c(&quot;black&quot;,&quot;red&quot;,&quot;green&quot;),xpd = T, cex = 2, y.intersp = 0.5)
</pre>
<p><img sizes="(max-width: 1000px) 100vw, 1000px" data-image-description="" data-orig-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot02.png?w=456&amp;ssl=1" data-orig-size="1000,800" data-image-title="rplot02" data-comments-opened="1" class="wp-image-1274" data-recalc-dims="1" data-attachment-id="1274" data-permalink="https://poissonisfish.wordpress.com/2017/01/23/principal-component-analysis-in-r/rplot02-2/" alt="rplot02" data-medium-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot02.png?w=456&amp;ssl=1?w=300" data-large-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot02.png?w=456&amp;ssl=1?w=1000" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot02.png?w=456&ssl=1"></p>
<p>Among other things, we observe correlations between variables (<em>e.g.</em> total phenols and flavonoids), and occasionally the two-dimensional separation of the three cultivars (<em>e.g.</em> using alcohol % and the OD ratio).</p>
<p>If its hard enough looking into all pairwise interactions in a set of 13 variables, let alone in sets of hundreds or thousands of variables. In these instances PCA is of great help. Let’s give it a try in this data set:</p>
<pre>
dev.off() # clear the format from the previous plot
winePCA &lt;- prcomp(scale(wine[,-1]))
plot(winePCA$x[,1:2], col = wineClasses)
</pre>
<p><img sizes="(max-width: 400px) 100vw, 400px" data-image-description="" data-orig-file="https://i2.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot.png?w=456&amp;ssl=1" data-orig-size="400,400" data-image-title="rplot" data-comments-opened="1" class="wp-image-1123" data-recalc-dims="1" data-attachment-id="1123" data-permalink="https://poissonisfish.wordpress.com/2017/01/23/principal-component-analysis-in-r/rplot-2/" alt="rplot" data-medium-file="https://i2.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot.png?w=456&amp;ssl=1?w=300" data-large-file="https://i2.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot.png?w=456&amp;ssl=1?w=400" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i2.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot.png?w=456&ssl=1"></p>
<p>Three lines of code and we see a clear separation among grape vine cultivars. In addition, the data points are evenly scattered over relatively narrow ranges in both PCs. We could next investigate which parameters contribute the most to this separation and how much variance is explained by each PC, but I will leave it for <strong>pcaMethods</strong>. We will now repeat the procedure after introducing an outlier in place of the 10th observation.</p>
<pre>
wineOutlier &lt;- wine
wineOutlier[10,] &lt;- wineOutlier[10,]*10 # change the 10th obs. into an extreme one by multiplying its profile by 10
outlierPCA &lt;- prcomp(scale(wineOutlier[,-1]))
plot(outlierPCA$x[,1:2], col = wineClasses)
</pre>
<p><img sizes="(max-width: 400px) 100vw, 400px" data-image-description="" data-orig-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot011.png?w=456&amp;ssl=1" data-orig-size="400,400" data-image-title="rplot01" data-comments-opened="1" class="wp-image-1138" data-recalc-dims="1" data-attachment-id="1138" data-permalink="https://poissonisfish.wordpress.com/2017/01/23/principal-component-analysis-in-r/rplot01/" alt="rplot01" data-medium-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot011.png?w=456&amp;ssl=1?w=300" data-large-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot011.png?w=456&amp;ssl=1?w=400" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot011.png?w=456&ssl=1"></p>
<p>As expected, the huge variance stemming from the separation of the 10th observation from the core of all other samples is fully absorbed by the first PC. The outlying sample becomes plain evident.</p>
<h3></h3>
<h3>PCA of the wine data set with pcaMethods</h3>
<p>We will now turn to <strong>pcaMethods</strong>, a compact suite of PCA tools. First you will need to install it from the Bioconductor:</p>
<pre>
source(&quot;https://bioconductor.org/biocLite.R&quot;)
biocLite(&quot;pcaMethods&quot;)
library(pcaMethods)
</pre>
<p> </p>
<p>There are three mains reasons why I use <strong>pcaMethods</strong> so extensively:</p>
<ul>
<li>Besides SVD, it provides several different methods (bayesian PCA, probabilistic PCA, robust PCA, to name a few)</li>
<li>Some of these algorithms tolerate and impute missing values</li>
<li>The object structure and plotting capabilities are user-friendly</li>
</ul>
<p>All information available about the package can be found <a target="_blank" rel="nofollow" href="http://bioconductor.org/packages/release/bioc/html/pcaMethods.html">here</a>. I will now simply show the joint scores-loadings plots, but still encourage you to explore it further.</p>
<p>I will select the default SVD method to reproduce our previous PCA result, with the same scaling strategy as before (UV, or unit-variance, as executed by <strong>scale</strong>). The argument <strong>scoresLoadings</strong> gives you control over printing scores, loadings, or both jointly as right next. The standard graphical parameters (<em>e.g.</em> <strong>cex</strong>, <strong>pch</strong>, <strong>col</strong>) preceded by either letters <strong>s</strong> or <strong>l</strong> control the aesthetics in the scores or loadings plots, respectively.</p>
<pre>
winePCAmethods &lt;- pca(wine[,-1], scale = &quot;uv&quot;, center = T, nPcs = 2, method = &quot;svd&quot;)
slplot(winePCAmethods, scoresLoadings = c(T,T), scol = wineClasses)
</pre>
<p><img sizes="(max-width: 799px) 100vw, 799px" data-image-description="" data-orig-file="https://i1.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot041.png?w=456&amp;ssl=1" data-orig-size="799,403" data-image-title="rplot04" data-comments-opened="1" class="wp-image-1388" data-recalc-dims="1" data-attachment-id="1388" data-permalink="https://poissonisfish.wordpress.com/2017/01/23/principal-component-analysis-in-r/rplot04/" alt="rplot04" data-medium-file="https://i1.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot041.png?w=456&amp;ssl=1?w=300" data-large-file="https://i1.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot041.png?w=456&amp;ssl=1?w=799" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i1.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot041.png?w=456&ssl=1"></p>
<p>So firstly, we have a faithful reproduction of the previous PCA plot. Then, having the loadings panel on its right side, we can claim that</p>
<ul>
<li>Wine from Cv2 (red) has a lighter color intensity, lower alcohol %, a greater OD ratio and hue, compared to the wine from Cv1 and Cv3.</li>
<li>Wine from Cv3 (green) has a higher content of malic acid and non-flavanoid phenols, and a higher alkalinity of ash compared to the wine from Cv1 (black)</li>
</ul>
<p>Finally, although the variance jointly explained by the first two PCs is printed by default (55.41%), it might be more informative consulting the variance explained in individual PCs. We can call the structure of <strong>winePCAmethods</strong>, inspect the slots and print those of interest, since there is a lot of information contained. The variance explained per component is stored in a slot named <strong>R2</strong>.</p>
<pre>
str(winePCAmethods) # slots are marked with @
winePCAmethods@R2
</pre>
<p>Seemingly, PC1 and PC2 explain 36.2% and 19.2% of the variance in the wine data set, respectively. PCAs of data exhibiting strong effects (such as the outlier example given above) will likely result in the sequence of PCs showing an abrupt drop in the variance explained. <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/stats/html/screeplot.html">Screeplots</a> are helpful in that matter, and allow you determining how much variance you can put into a principal component regression (PCR), for example, which is exactly what we will try next.</p>
<h3></h3>
<h3>PCR with the housing data set</h3>
<p>Now we will tackle a regression problem using PCR. I will use an old housing data set also deposited in the UCI MLR. Again according to its <a target="_blank" rel="nofollow" href="http://archive.ics.uci.edu/ml/datasets/Housing">documentation</a>, these data consist of 14 variables and 504 records from distinct towns somewhere in the US. To perform PCR all we need is conduct PCA and feed the scores of <img title="n " alt="n " class="latex" src="https://s0.wp.com/latex.php?latex=n+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> PCs to a OLS. Let’s try predicting the median value of owner-occupied houses in thousands of dollars (<strong>MEDV</strong>) using the first three PCs from a PCA.</p>
<pre>
houses &lt;- read.table(&quot;http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data&quot;,header = F, na.string = &quot;?&quot;)
colnames(houses) &lt;- c(&quot;CRIM&quot;, &quot;ZN&quot;, &quot;INDUS&quot;,&quot;CHAS&quot;,&quot;NOX&quot;,&quot;RM&quot;,&quot;AGE&quot;,&quot;DIS&quot;,&quot;RAD&quot;,&quot;TAX&quot;,&quot;PTRATIO&quot;,&quot;B&quot;,&quot;LSTAT&quot;,&quot;MEDV&quot;)

# Perform PCA
pcaHouses &lt;- prcomp(scale(houses[,-14]))
scoresHouses &lt;- pcaHouses$x

# Fit lm using the first 3 PCs
modHouses &lt;- lm(houses$MEDV ~ scoresHouses[,1:3])
summary(modHouses)
</pre>
<p>The printed summary shows two important pieces of information. Firstly, the three estimated coefficients (plus the intercept) are considered significant (<img title="H_0 : \beta = 0 " alt="H_0 : \beta = 0 " class="latex" src="https://s0.wp.com/latex.php?latex=H_0+%3A+%5Cbeta+%3D+0+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">). Second, the predictability as defined by the <img title="R^2 " alt="R^2 " class="latex" src="https://s0.wp.com/latex.php?latex=R%5E2+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> (coefficient of determination, in most cases the same as the squared Pearson correlation coefficient) was 0.63.</p>
<p>Next we will compare this simple model to a OLS model featuring all 14 variables, and finally compare the observed vs. predicted MEDV plots from both models. Note that in the <strong>lm</strong> syntax, the response <img title="Y " alt="Y " class="latex" src="https://s0.wp.com/latex.php?latex=Y+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> is given to the left of the tilde and the set of predictors <img title="X " alt="X " class="latex" src="https://s0.wp.com/latex.php?latex=X+&bg=f9f9f9&%23038;fg=555555&%23038;s=1"> to the right. Moreover, provided there is an argument for <strong>data</strong> you can circumvent the need for typing all variable names for a full model (<img title="x_1 + x_2 + x_3 + ... x_p " alt="x_1 + x_2 + x_3 + ... x_p " class="latex" src="https://s0.wp.com/latex.php?latex=x_1+%2B+x_2+%2B+x_3+%2B+...+x_p+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">), and simply use <strong>. </strong>.</p>
<pre>
# Fit lm using all 14 vars
modHousesFull &lt;- lm(MEDV ~ ., data = houses)
summary(modHousesFull) # R2 = 0.741

# Compare obs. vs. pred. plots
par(mfrow = c(1,2))
plot(houses$MEDV, predict(modHouses), xlab = &quot;Observed MEDV&quot;, ylab = &quot;Predicted MEDV&quot;, main = &quot;PCR&quot;, abline(a = 0, b = 1, col = &quot;red&quot;))
plot(houses$MEDV, predict(modHousesFull), xlab = &quot;Observed MEDV&quot;, ylab = &quot;Predicted MEDV&quot;, main = &quot;Full model&quot;, abline(a = 0, b = 1, col = &quot;red&quot;))
</pre>
<p><img sizes="(max-width: 800px) 100vw, 800px" data-image-description="" data-orig-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot03.png?w=456&amp;ssl=1" data-orig-size="800,400" data-image-title="rplot03" data-comments-opened="1" class="wp-image-1316" data-recalc-dims="1" data-attachment-id="1316" data-permalink="https://poissonisfish.wordpress.com/2017/01/23/principal-component-analysis-in-r/rplot03-2/" alt="rplot03" data-medium-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot03.png?w=456&amp;ssl=1?w=300" data-large-file="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot03.png?w=456&amp;ssl=1?w=800" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i0.wp.com/poissonisfish.files.wordpress.com/2017/01/rplot03.png?w=456&ssl=1"></p>
<p>Here the full model displays a slight improvement in fit (<img title="R^2 = 0.74 " alt="R^2 = 0.74 " class="latex" src="https://s0.wp.com/latex.php?latex=R%5E2+%3D+0.74+&bg=f9f9f9&%23038;fg=555555&%23038;s=1">). The high significance of most coefficient estimates is suggestive of a well-designed experiment. Nevertheless, it is notable that such a reduction of 13 down to three covariates still yields an accurate model.</p>
<p>Just as a side note, you probably noticed both models underestimated the MEDV in towns with MEVD worth 50,000 dollars. My guess is that missing values were set to MEVD = 50.</p>
<p>Concluding,</p>
<ul>
<li>The SVD algorithm is founded on fundamental properties of linear algebra including matrix diagonalization. SVD-based PCA takes part of its solution and retains a reduced number of orthogonal covariates that explain as much variance as possible.</li>
<li>Use PCA when handling high-dimensional data. It is insensitive to correlation among variables and efficient in detecting sample outliers.</li>
<li>If you plan to use PCA results for subsequent analyses all care should be undertaken in the process. Although typically outperformed by numerous methods, PCR still benefits from interpretability and can be effective in many settings.</li>
</ul>
<p>All feedback from these tutorials is very welcome, please enter the Contact tab and leave your comments. I do also appreciate suggestions. Enjoy!</p>
<p> </p>
<p> </p>
<p> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/poissonisfish.wordpress.com/572/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/poissonisfish.wordpress.com/572/" border="0"></a> <img alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" border="0" data-recalc-dims="1"></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="https://poissonisfish.wordpress.com/2017/01/23/principal-component-analysis-in-r/"> poissonisfish</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/7cGsv1XqhCI"></a></div>