<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/obtaining-the-number-of-components-from-cross-validation-of-principal-components-regression/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://statisticaloddsandends.wordpress.com/2018/10/15/obtaining-the-number-of-components-from-cross-validation-of-principal-components-regression/"> R – Statistical Odds &amp; Ends</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><a rel="nofollow" href="https://en.wikipedia.org/wiki/Principal_component_regression" target="_blank">Principal components (PC) regression</a> is a common dimensionality reduction technique in supervised learning. The R lab for PC regression in James et al.’s <a rel="nofollow" href="https://www-bcf.usc.edu/~gareth/ISL/ISLR%20Seventh%20Printing.pdf" target="_blank"><em>Introduction to Statistical Learning</em></a> is a popular intro for how to perform PC regression in R: it is on p256-257 of the book (p270-271 of the PDF).</p>
<p>As in the lab, the code below runs PC regression on the <code>Hitters</code> data to predict <code>Salary</code>:</p>
<pre> 
library(ISLR)
data(Hitters)
library(pls)
set.seed(2)
pcr.fit = pcr(Salary ~ ., data = Hitters, scale = TRUE, 
    validation = &quot;CV&quot;) 
</pre>
<p>The result of the fit can be seen by calling <code>summary(pcr.fit)</code>:</p>
<div data-shortcode="caption" class="wp-caption"><img data-orig-file="https://statisticaloddsandends.files.wordpress.com/2018/09/summary-pcr.png?w=456" data-orig-size="520,290" data-medium-file="https://statisticaloddsandends.files.wordpress.com/2018/09/summary-pcr.png?w=456?w=300" data-attachment-id="442" data-image-meta="{" data-image-description="" src="https://statisticaloddsandends.files.wordpress.com/2018/09/summary-pcr.png?w=456" data-permalink="https://statisticaloddsandends.wordpress.com/2018/10/15/obtaining-the-number-of-components-from-cross-validation-of-principal-components-regression/summary-pcr/" data-large-file="https://statisticaloddsandends.files.wordpress.com/2018/09/summary-pcr.png?w=456?w=520" alt="" data-image-title="summary pcr" sizes="(max-width: 520px) 100vw, 520px" data-comments-opened="1" class="wp-image-442" data-recalc-dims="1">
<p class="wp-caption-text">Taken from ISLR Section 6.7.1.</p>
</div>
<p>When predicting the response on test data, we would choose the number of components based on which model gave the smallest cross validation (CV) error. What is strange in this lab is that they do not pull out the “optimal” number of components via code! Instead, they make a plot of CV error, visually identified the number of components giving smallest CV error (7 in this case), <strong>then predicted on the test set with 7 hardcoded</strong>. Here is their code:</p>
<pre> 
# get indices for training set
set.seed(1)
train = sample(1:nrow(Hitters), nrow(Hitters) / 2)

# perform PC regression with CV on just the training set
set.seed(1)
pcr.fit = pcr(Salary ~ ., data = Hitters, subset = train, 
              scale = TRUE, validation = &quot;CV&quot;)

# plot of CV error
validationplot(pcr.fit, val.type = &quot;MSEP&quot;)
# visual inspection of plot: 7 gives min CV error!

# predict on test data
pcr.pred = predict(pcr.fit, x[-train,], ncomp = 7)
</pre>
<p>This method of identifying the number of components with smallest CV error is not ideal. Here is the output of <code>validationplot</code>:</p>
<p><img data-orig-file="https://statisticaloddsandends.files.wordpress.com/2018/10/validationplot.png?w=456" data-orig-size="358,336" data-medium-file="https://statisticaloddsandends.files.wordpress.com/2018/10/validationplot.png?w=456?w=300" data-attachment-id="481" data-image-meta="{" data-image-description="" src="https://statisticaloddsandends.files.wordpress.com/2018/10/validationplot.png?w=456" data-permalink="https://statisticaloddsandends.wordpress.com/2018/10/15/obtaining-the-number-of-components-from-cross-validation-of-principal-components-regression/validationplot/" data-large-file="https://statisticaloddsandends.files.wordpress.com/2018/10/validationplot.png?w=456?w=358" alt="" data-image-title="validationplot" sizes="(max-width: 358px) 100vw, 358px" data-comments-opened="1" class="wp-image-481" data-recalc-dims="1"></p>
<p>I think it’s actually quite hard to tell by eye that 7 principal components is the best choice! In any case, we might want a programmatic way of finding the best number of components (according to CV) so that we don’t have to stop our data modeling pipeline to make a visual check.</p>
<p>It turns out that we can use the <code>RMSEP</code> function to pull out the CV and adjCV table that we saw in <code>summary(pcr.fit)</code>:</p>
<pre>
RMSEP(pcr.fit)
#        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps  9 comps  10 comps
# CV             488    386.2    387.3    387.0    387.9    390.9    383.8    382.5    388.0    388.1     385.2
# adjCV          488    385.1    386.2    385.7    386.5    389.5    381.1    380.5    385.7    385.7     382.7
#        11 comps  12 comps  13 comps  14 comps  15 comps  16 comps  17 comps  18 comps  19 comps
# CV        385.9     387.1     392.7     408.1     395.3     386.6     394.0     399.6     405.6
# adjCV     383.3     384.5     390.1     404.6     391.3     382.8     389.8     394.8     400.5
</pre>
<p>To get at the individual values in the table, we need to access <code>val</code> key and subset appropriately:</p>
<pre>
cverr &lt;- RMSEP(pcr.fit)$val[1,,]
imin &lt;- which.min(cverr) - 1
</pre>
<p>Note that we have to subtract 1, since the model with no principal components is included in the output of <code>RMSEP(pcr.fit)</code> as well. If we wish to get the best number of PCs based on adjusted CV instead, we could amend that first line of code to <code>RMSEP(pcr.fit)$val[1,,]</code>.</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://statisticaloddsandends.wordpress.com/2018/10/15/obtaining-the-number-of-components-from-cross-validation-of-principal-components-regression/"> R – Statistical Odds &amp; Ends</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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/5M_-twhWEjU" height="1"></a></div>