<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/principal-components-regression-pt-3-picking-the-number-of-components/"></div></div><div><a data-url="http://www.r-bloggers.com/principal-components-regression-pt-3-picking-the-number-of-components/" data-counturl="http://www.r-bloggers.com/principal-components-regression-pt-3-picking-the-number-of-components/" data-count="horizontal" data-text="Principal Components Regression, Pt. 3: Picking the Number of Components" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.win-vector.com/blog/2016/05/pcr_part3_pickk/"> R – Win-Vector Blog</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>In our <a target="_blank" href="http://www.win-vector.com/blog/2016/05/pcr_part2_yaware">previous note</a> we demonstrated <em>Y</em>-Aware PCA and other <em>y</em>-aware approaches to dimensionality reduction in a predictive modeling context, specifically Principal Components Regression (PCR). For our examples, we selected the appropriate number of principal components by eye. In this note, we will look at ways to select the appropriate number of principal components in a more automated fashion.</p>
<p><span></span></p>
<p>Before starting the discussion, let’s quickly redo our <em>y</em>-aware PCA. Please refer to our <a target="_blank" href="http://www.win-vector.com/blog/2016/05/pcr_part2_yaware">previous post</a> for a full discussion of this data set and this approach.</p>
<pre><code>#
# make data
#
set.seed(23525)
dTrain &lt;- mkData(1000)
dTest &lt;- mkData(1000)

#
# design treatment plan
#
treatmentsN &lt;- designTreatmentsN(dTrain,
                                 setdiff(colnames(dTrain),'y'),'y',
                                 verbose=FALSE)

#
# prepare the treated frames, with y-aware scaling
#
examplePruneSig = 1.0 
dTrainNTreatedYScaled &lt;- prepare(treatmentsN,dTrain,
                                 pruneSig=examplePruneSig,scale=TRUE)
dTestNTreatedYScaled &lt;- prepare(treatmentsN,dTest,
                                pruneSig=examplePruneSig,scale=TRUE)

#
# do the principal components analysis
#
vars &lt;- setdiff(colnames(dTrainNTreatedYScaled),'y')
# prcomp defaults to scale. = FALSE, but we already 
# scaled/centered in vtreat- which we don't want to lose.
dmTrain &lt;- as.matrix(dTrainNTreatedYScaled[,vars])
dmTest &lt;- as.matrix(dTestNTreatedYScaled[,vars])
princ &lt;- prcomp(dmTrain, center = FALSE, scale. = FALSE)</code></pre>
<p>If we examine the magnitudes of the resulting singular values, we see that we should use from two to five principal components for our analysis. In fact, as we showed in the previous post, the first two singular values accurately capture the two unobservable processes that contribute to <em>y</em>, and a linear model fit to these two components captures most of the explainable variance in the data, both on training and on hold-out data.</p>
<p><img title="" alt="" src="http://i2.wp.com/www.win-vector.com/dfiles/YAwarePCR_pickK_files/figure-html/singularvalues-1.png?w=660" data-recalc-dims="1"></p>
<p>We picked the number of principal components to use by eye; but it’s tricky to implement code based on the strategy &quot;look for a knee in the curve.&quot; So how might we automate picking the appropriate number of components in a reliable way?</p>
<div>
<h2><em>X</em>-Only Approaches</h2>
<p>Jackson (1993) and Peres-Neto, et.al. (2005) are two excellent surveys and evaluations of the different published approaches to picking the number of components in standard PCA. Those methods include:</p>
<ol>
<li>Look for a &quot;knee in the curve&quot; — the approach we have taken, visually.</li>
<li><em>Only</em> for data that has been scaled to unit variance: keep the components corresponding to singular values greater than 1.</li>
<li>Select enough components to cover some fixed fraction (generally 95%) of the observed variance. This is the approach taken by <code>caret::preProcess.</code></li>
<li>Perform a statistical test to see which singular values are larger than we would expect from an appropriate null hypothesis or noise process.</li>
</ol>
<p>The papers also cover other approaches, as well as different variations of the above.</p>
<p>Kabakoff (<em>R In Action</em>, 2nd Edition, 2015) suggests comparing the magnitudes of the singular values to those extracted from random matrices of the same shape as the original data. Let’s assume that the original data has <em>k</em> variables, and that PCA on the original data extracts the <em>k</em> singular values s<sub>i</sub> and the <em>k</em> principal components PC<sub>i</sub>.To pick the appropriate number of principal components:</p>
<ol>
<li>For a chosen number of iterations, <em>N</em> (choose <em>N</em> >> <em>k</em>):</li>
</ol>
<ul>
<li>Generate a random matrix of the correct size</li>
<li>Do PCA and extract the singular values</li>
</ul>
<ol start="2">
<li>Then for each of the <em>k</em> principal components:</li>
</ol>
<ul>
<li>Find the mean of the <em>i</em>th singular value, r<sub>i</sub></li>
<li>If s<sub>i</sub> > r<sub>i</sub>, then keep PC<sub>i</sub></li>
</ul>
<p>The idea is that if there is more variation in a given direction than you would expect at random, then that direction is probably meaningful. If you assume that higher variance directions are more useful than lower variance directions (the usual assumption), then one handy variation is to find the first <em>i</em> such that s<sub>i</sub> &lt; r<sub>i</sub>, and keep the first <em>i</em>-1 principal components.</p>
<p>This approach is similar to what the authors of the survey papers cited above refer to as the <em>broken-stick method</em>. In their research, the broken-stick method was among the best performing approaches for a variety of simulated and real-world examples.</p>
<p>With the proper adjustment, all of the above heuristics work as well in the <em>y</em>-adjusted case as they do with traditional <em>x</em>-only PCA.</p>
</div>
<div>
<h2>A <em>Y</em>-Aware Approach: The Permutation Test</h2>
<p>Since in our case we know <em>y</em>, we can — and should — take advantage of this information. We will use a variation of the broken-stick method, but rather than comparing our data to a random matrix, we will compare our data to alternative datasets where <em>x</em> has no relation to <em>y</em>. We can do this by randomly permuting the <em>y</em> values. This preserves the structure of <em>x</em> — that is, the correlations and relationships of the <em>x</em> variables to each other — but it changes the units of the problem, that is, the <em>y</em>-aware scaling. We are testing whether or not a given principal component appears more meaningful in a metric space induced by the true <em>y</em> than it does in a random metric space, one that preserves the distribution of <em>y</em>, but not the relationship of <em>y</em> to <em>x</em>.</p>
<p>You can read a more complete discussion of permutation tests and their application to variable selection (significance pruning) in <a target="_blank" href="http://www.win-vector.com/blog/2015/08/how-do-you-know-if-your-data-has-signal/">this post</a>.</p>
<p>In our example, we’ll use <em>N</em>=100, and rather than using the means of the singular values from our experiments as the thresholds, we’ll use the 98th percentiles. This represents a threshold value that is likely to be exceeded by a singular value induced in a random space only 1/(the number of variables) (1/50=0.02) fraction of the time.</p>
<pre><code>#
# Resample y, do y-aware PCA, 
# and return the singular values
#
getResampledSV = function(data,yindices) {
  # resample y
  data$y = data$y[yindices]
  
  # treatment plan
  treatplan = vtreat::designTreatmentsN(data, 
                                setdiff(colnames(data), 'y'), 
                                'y', verbose=FALSE)
  # y-aware scaling
  dataTreat = vtreat::prepare(treatplan, data, pruneSig=1, scale=TRUE)
  
  # PCA
  vars = setdiff(colnames(dataTreat), 'y')
  dmat = as.matrix(dataTreat[,vars])
  princ = prcomp(dmat, center=FALSE, scale=FALSE)
  
  # return the magnitudes of the singular values
  princ$sdev
}

#
# Permute y, do y-aware PCA, 
# and return the singular values
#
getPermutedSV = function(data) {
  n = nrow(data)
  getResampledSV(data,sample(n,n,replace=FALSE))
}

#
# Run the permutation tests and collect the outcomes
#
niter = 100 # should be >> nvars
nvars = ncol(dTrain)-1
# matrix: 1 column for each iter, nvars rows
svmat = vapply(1:niter, FUN=function(i) {getPermutedSV(dTrain)}, numeric(nvars))
rownames(svmat) = colnames(princ$rotation) # rows are principal components
colnames(svmat) = paste0('rep',1:niter) # each col is an iteration

# plot the distribution of values for the first singular value
# compare it to the actual first singular value
ggplot(as.data.frame(t(svmat)), aes(x=PC1)) + 
  geom_density() + geom_vline(xintercept=princ$sdev[[1]], color=&quot;red&quot;) +
  ggtitle(&quot;Distribution of magnitudes of first singular value, permuted data&quot;)</code></pre>
<p><img title="" alt="" src="http://i2.wp.com/www.win-vector.com/dfiles/YAwarePCR_pickK_files/figure-html/permutation-1.png?w=660" data-recalc-dims="1"></p>
<p>Here we show the distribution of the magnitude of the first singular value on the permuted data, and compare it to the magnitude of the actual first singular value (the red vertical line). We see that the actual first singular value is far larger than the magnitude you would expect from data where <em>x</em> is not related to <em>y</em>. Let’s compare all the singular values to their permutation test thresholds. The dashed line is the mean value of each singular value from the permutation tests; the shaded area represents the 98th percentile.</p>
<pre><code># transpose svmat so we get one column for every principal component
# Get the mean and empirical confidence level of every singular value
as.data.frame(t(svmat)) %>% dplyr::summarize_each(funs(mean)) %>% as.numeric() -> pmean
confF &lt;- function(x) as.numeric(quantile(x,1-1/nvars))
as.data.frame(t(svmat)) %>% dplyr::summarize_each(funs(confF)) %>% as.numeric() -> pupper

pdata = data.frame(pc=seq_len(length(pmean)), magnitude=pmean, upper=pupper)

# we will use the first place where the singular value falls 
# below its threshold as the cutoff.
# Obviously there are multiple comparison issues on such a stopping rule,
# but for this example the signal is so strong we can ignore them.
below = which(princ$sdev &lt; pdata$upper)
lastSV = below[[1]] - 1</code></pre>
<p><img title="" alt="" src="http://i1.wp.com/www.win-vector.com/dfiles/YAwarePCR_pickK_files/figure-html/replot-1.png?w=660" data-recalc-dims="1"></p>
<p>This test suggests that we should use 5 principal components, which is consistent with what our eye sees. This is perhaps not the &quot;correct&quot; knee in the graph, but it is undoubtably a knee.</p>
<div>
<h3>Bootstrapping</h3>
<p>Empirically estimating the quantiles from the permuted data so that we can threshold the non-informative singular values will have some undesirable bias and variance, especially if we do not perform enough experiment replications. This suggests that instead of estimating quantiles ad-hoc, we should use a systematic method: The Bootstrap. Bootstrap replication breaks the input to output association by re-sampling with replacement rather than using permutation, but comes with built-in methods to estimate bias-adjusted confidence intervals. The methods are fairly technical, and on this dataset the results are similar, so we don’t show them here, although the code is available in the <a target="_blank" href="https://github.com/WinVector/Examples/blob/master/PCR/YAwarePCR_pickK.Rmd">R markdown document used to produce this note</a>.</p>
</div>
</div>
<div>
<h2>Significance Pruning</h2>
<p>Alternatively, we can treat the principal components that we extracted via <em>y</em>-aware PCA simply as transformed variables — which is what they are — and significance prune them in the standard way. As our <a target="_blank" href="http://www.win-vector.com/blog/2015/08/how-do-you-know-if-your-data-has-signal/">article on significance pruning</a> discusses, we can estimate the significance of a variable by fitting a one variable model (in this case, a linear regression) and looking at that model’s significance value. You can pick the pruning threshold by considering the rate of false positives that you are willing to tolerate; as a rule of thumb, we suggest one over the number of variables.</p>
<p>In regular significance pruning, you would take any variable with estimated significance value lower than the threshold. Since in the PCR situation we presume that the variables are ordered from most to least useful, you can again look for the first position <em>i</em> where the variable appears insignificant, and use the first <em>i</em>-1 variables.</p>
<p>We’ll use <code>vtreat</code> to get the significance estimates for the principal components. We’ll use one over the number of variables (1/50 = 0.02) as the pruning threshold.</p>
<pre><code># get all the principal components
# not really a projection as we took all components!
projectedTrain &lt;- as.data.frame(predict(princ,dTrainNTreatedYScaled),
                                 stringsAsFactors = FALSE)
vars = colnames(projectedTrain)
projectedTrain$y = dTrainNTreatedYScaled$y

# designing the treatment plan for the transformed data
# produces a data frame of estimated significances
tplan = designTreatmentsN(projectedTrain, vars, 'y', verbose=FALSE)

threshold = 1/length(vars)
scoreFrame = tplan$scoreFrame
scoreFrame$accept = scoreFrame$sig &lt; threshold

# pick the number of variables in the standard way:
# the number of variables that pass the significance prune
nPC = sum(scoreFrame$accept)</code></pre>
<p><img title="" alt="" src="http://i0.wp.com/www.win-vector.com/dfiles/YAwarePCR_pickK_files/figure-html/plotsig-1.png?w=660" data-recalc-dims="1"></p>
<p>Significance pruning picks 2 principal components, again consistent with our visual assessment. This time, we picked the correct knee: as we saw in the previous post, the first two principal components were sufficient to describe the explainable structure of the problem.</p>
</div>
<div>
<h2>Conclusion</h2>
<p>Since one of the purposes of PCR/PCA is to discover the underlying structure in the data, it’s generally useful to examine the singular values and the variable loadings on the principal components. However an analysis should also be repeatable, and hence, automatable, and it’s not straightforward to automate something as vague as &quot;look for a knee in the curve&quot; when selecting the number of principal components to use. We’ve covered two ways to programatically select the appropriate number of principal components in a predictive modeling context.</p>
<p>To conclude this entire series, here is our recommended best practice for principal components regression:</p>
<ol>
<li>Significance prune the candidate input variables.</li>
<li>Perform a Y-Aware principal components analysis.</li>
<li>Significance prune the resulting principal components.</li>
<li>Regress.</li>
</ol>
<p>Thanks to <a target="_blank" href="https://plus.google.com/+CyrilPernet">Cyril Pernet</a>, who blogs at <a target="_blank" href="http://neurostatscyrilpernet.blogspot.co.uk/"><em>NeuroImaging and Statistics</em></a>, for requesting this follow-up post and pointing us to the Jackson reference.</p>
<div>
<h3>References</h3>
<ul>
<li>
<p>Jackson, Donald A. &quot;Stopping Rules in Principal Components Analysis: A Comparison of Heuristical and Statistical Approaches&quot;, <em>Ecology</em> Vol 74, no. 8, 1993.</p>
</li>
<li>
<p>Kabacoff, Robert I. <em>R In Action</em>, 2nd edition, Manning, 2015.</p>
</li>
<li>
<p>Efron, Bradley and Robert J. Tibshirani. <em>An Introduction to the Bootstrap</em>, Chapman and Hall/CRC, 1998.</p>
</li>
<li>
<p>Peres-Neto, Pedro, Donald A. Jackson and Keith M. Somers. &quot;How many principal components? Stopping rules for determining the number of non-trivial axes revisited&quot;, <em>Computational Statistics &amp; Data Analysis</em>, Vol 49, no. 4, 2005.</p>
</li>
</ul>
</div>
</div>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/principal-components-regression-pt-3-picking-the-number-of-components/"></div></div><div><a data-url="http://www.r-bloggers.com/principal-components-regression-pt-3-picking-the-number-of-components/" data-counturl="http://www.r-bloggers.com/principal-components-regression-pt-3-picking-the-number-of-components/" data-count="vertical" data-text="Principal Components Regression, Pt. 3: Picking the Number of Components" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<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.win-vector.com/blog/2016/05/pcr_part3_pickk/"> R – Win-Vector Blog</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://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="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/f2UK4JDXEoo"></a></div>