 |
| See Video ⮞ ☝ |
AGRON Stats
May 26, 2020
Introduction
Principal component analysis (PCA) reduces a set of correlated numeric variables to a smaller set of orthogonal components. Each component is a linear combination of the original variables, ordered by the amount of variance it captures. PCA is exploratory and its usefulness depends on the variables, preprocessing, and scientific question.
PCA will generate as many components as there are variables. When the first few components retain a substantial share of variation, they can provide a compact summary for visualization or downstream exploration. The retained components do not automatically preserve every feature relevant to a later prediction or causal analysis.
PCA is often used as a first step for further multivariate data analysis procedures like:
- Cluster analysis
- Multiple regression
- Discriminant analysisPCA can be computed from a covariance or correlation matrix. The choice determines how variable scales contribute: standardizing variables gives them equal variance before the decomposition, while using the covariance matrix preserves their original scale.
This example uses measurements from 200 banknotes, with the first 100 genuine and the remaining 100 counterfeit. Verify the group labels and row order in the source data rather than assuming observations are arranged as described.
Clear objects from the workspace only if you have saved anything needed for later. Closing graphics devices and clearing the console are optional and do not affect the PCA.
Import data set
Import the CSV file with read.csv(), providing the correct path and confirming that the first row contains variable names when header = TRUE. Inspect the column names and dimensions after import.
Use head() and str() to verify the data and identify the categorical group variable. PCA must be calculated on numeric measurements, not on identifiers or group labels.
# note.type Length Left Right Bottom Top Diagonal
# 1 counterfeit 214.8 131.0 131.1 9.0 9.7 141.0
# 2 counterfeit 214.6 129.7 129.7 8.1 9.5 141.7
# 3 counterfeit 214.8 129.7 129.7 8.7 9.6 142.2
# 4 counterfeit 214.8 129.7 129.6 7.5 10.4 142.0
# 5 counterfeit 215.0 129.6 129.7 10.4 7.7 141.8
# 6 counterfeit 215.7 130.8 130.5 9.0 10.1 141.4Standardize the data set
Centering and scaling should be chosen deliberately. Variables with larger units or variances can dominate a covariance-based PCA. Standardization prevents this scale dominance, but it also gives low-variance variables equal weight and may amplify measurement noise.
The scale() function can center columns and divide them by their standard deviations. Exclude the categorical identifier with an explicit column selection, and check that no selected column is constant or otherwise unsuitable for scaling.
Inspect the transformed values to confirm that centering and scaling were applied as intended. Standardized columns should have means near zero and, where defined, standard deviations near one.
# Length Left Right Bottom Top Diagonal
# [1,] -0.2549435 2.433346 2.8299417 -0.2890067 -1.1837648 0.4482473
# [2,] -0.7860757 -1.167507 -0.6347880 -0.9120152 -1.4328473 1.0557460
# [3,] -0.2549435 -1.167507 -0.6347880 -0.4966762 -1.3083061 1.4896737
# [4,] -0.2549435 -1.167507 -0.8822687 -1.3273542 -0.3119759 1.3161027
# [5,] 0.2761888 -1.444496 -0.6347880 0.6801176 -3.6745902 1.1425316
# [6,] 2.1351516 1.879368 1.3450576 -0.2890067 -0.6855997 0.7953894Extraction of components (eigenvalues)
Eigen values & vectors
A PCA decomposition can be computed with prcomp(), princomp(), or an eigendecomposition of a covariance or correlation matrix. Use a consistent preprocessing choice when comparing their outputs.
See other tuotials where these methods have been used and explained —
PCA_1andPCA_2
The example uses eigen() to decompose a matrix of variable covariances or correlations. Confirm which matrix is supplied: the eigenvalues and loadings will differ depending on whether variables were standardized.
# eigen() decomposition
# $values
# [1] 3.00030487 0.93562052 0.24341371 0.19465874 0.08521185 0.03551468
#
# $vectors
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] -0.04377427 0.01070966 -0.3263165 -0.5616918 0.75257278 0.09809807
# [2,] 0.11216159 0.07144697 -0.2589614 -0.4554588 -0.34680082 -0.76651197
# [3,] 0.13919062 0.06628208 -0.3447327 -0.4153296 -0.53465173 0.63169678
# [4,] 0.76830499 -0.56307225 -0.2180222 0.1861082 0.09996771 -0.02221711
# [5,] 0.20176610 0.65928988 -0.5566857 0.4506985 0.10190229 -0.03485874
# [6,] -0.57890193 -0.48854255 -0.5917628 0.2584483 -0.08445895 -0.04567946Eigenvectors define the directions of the principal axes, and the corresponding eigenvalues quantify variance along those directions. Component signs can be reversed without changing the PCA solution.
An eigenvector is a set of coefficients defining a direction in the feature space. These coefficients, often called loadings under a stated convention, help describe how each original variable contributes to a component.
An eigenvalue measures the variance captured along its associated principal component. Larger eigenvalues indicate components that account for more variance in the analyzed matrix.
The first two components are often plotted for exploration, but the number to retain should reflect variance, interpretability, and the analytical objective.
Repeat the decomposition for the standardized data when comparing a correlation-based PCA with a covariance-based PCA. Keep the two preprocessing choices clearly labeled.
# eigen() decomposition
# $values
# [1] 2.9455582 1.2780838 0.8690326 0.4497687 0.2686769 0.1888799
#
# $vectors
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] -0.006987029 0.81549497 -0.01768066 0.5746173 -0.0587961 -0.03105698
# [2,] 0.467758161 0.34196711 0.10338286 -0.3949225 0.6394961 0.29774768
# [3,] 0.486678705 0.25245860 0.12347472 -0.4302783 -0.6140972 -0.34915294
# [4,] 0.406758327 -0.26622878 0.58353831 0.4036735 -0.2154756 0.46235361
# [5,] 0.367891118 -0.09148667 -0.78757147 0.1102267 -0.2198494 0.41896754
# [6,] -0.493458317 0.27394074 0.11387536 -0.3919305 -0.3401601 0.63179849Standardization changes the relative contribution of variables and therefore can change the eigenvalues, loadings, and scores. Compare the solutions in the context of the measurement scales and scientific question, not as interchangeable versions.
Proportion of variance
Calculate each component’s proportion of variance by dividing its eigenvalue by the sum of all eigenvalues. This proportion is meaningful only for the matrix used to fit the PCA.
The cumulative proportion is the sum of variance proportions through each component. In this example, the reported values summarize the standardized analysis; recompute them if the variables or preprocessing change.
e.scaled$values # Variances
e.scaled$values/6 # proportion of variances
cumsum(e.scaled$values/6) # Cumulative proportion# [1] 2.9455582 1.2780838 0.8690326 0.4497687 0.2686769 0.1888799
# [1] 0.49092637 0.21301396 0.14483876 0.07496145 0.04477948 0.03147998
# [1] 0.4909264 0.7039403 0.8487791 0.9237405 0.9685200 1.0000000Computation of principal component vectors
Component scores can be obtained by projecting the centered data onto the eigenvectors. Use the same centering and scaling parameters used to fit the PCA; mixing raw and standardized data produces incompatible scores.
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 1.7430272 1.64669605 1.4201973 -2.7479691 0.003293759 -0.60202200
# [2,] -2.2686248 -0.53744461 0.5313151 -0.6573558 -0.158171742 -0.45654268
# [3,] -2.2717009 -0.10740754 0.7156191 -0.3408384 -0.453880889 0.04532905
# [4,] -2.2778385 -0.08743490 -0.6041176 -0.3918255 -0.282913485 0.05543875
# [5,] -2.6255397 0.03909779 3.1883837 0.4240168 -0.277502895 -0.72026433
# [6,] 0.7565089 3.08101359 0.7845117 -0.5980322 0.192757017 0.10529393Computation of vectors from stats package
The prcomp() function performs the same analysis within base R. Supply the numeric feature matrix and choose center and scale. to match the intended covariance- or correlation-based analysis.
# PC1 PC2 PC3 PC4 PC5 PC6
# [1,] -1.7430272 -1.64669605 -1.4201973 -2.7479691 0.003293759 0.60202200
# [2,] 2.2686248 0.53744461 -0.5313151 -0.6573558 -0.158171742 0.45654268
# [3,] 2.2717009 0.10740754 -0.7156191 -0.3408384 -0.453880889 -0.04532905
# [4,] 2.2778385 0.08743490 0.6041176 -0.3918255 -0.282913485 -0.05543875
# [5,] 2.6255397 -0.03909779 -3.1883837 0.4240168 -0.277502895 0.72026433
# [6,] -0.7565089 -3.08101359 -0.7845117 -0.5980322 0.192757017 -0.10529393The PC values are the same as obtained using eigen() function however the signs are reversed. Use summary() to inspect standard deviations, proportions of variance, and cumulative proportions. Confirm that these values agree with the chosen preprocessing.
# Importance of components:
# PC1 PC2 PC3 PC4 PC5 PC6
# Standard deviation 1.7163 1.1305 0.9322 0.67065 0.51834 0.43460
# Proportion of Variance 0.4909 0.2130 0.1448 0.07496 0.04478 0.03148
# Cumulative Proportion 0.4909 0.7039 0.8488 0.92374 0.96852 1.00000Display a screeplot
A scree plot displays eigenvalues or component variance and can help identify where additional components contribute relatively little. Interpret it alongside cumulative variance and component meaning.
Use screeplot() to display the fitted PCA object. State the type of scale shown and check the component order and labels.
A bar plot of component variances is another way to show their relative contributions. Avoid implying that a visual elbow is a definitive statistical selection rule.
par(mfrow=c(1, 3))
plot(e.scaled$values, xlab = "Index", ylab = "Lambda",
main = "Scree plot",
cex.lab = 1.2, cex.axis = 1.2, cex.main = 1.8)
screeplot(x= pc, type = "line", main = "Scree plot", cex.main = 1.8)
plot(pc, main = "Variance bar plot", cex.main = 1.8)Correlations between original variables and component scores can help describe how variables align with the components. The calculation must use the same observations and preprocessing as the fitted PCA.
Most of the variation is due to the first two prinicple components (70%).Correlation of original variables with PCs
Center the numeric variables using the same means used in the fitted analysis. If the PCA was standardized, divide by the same standard deviations before calculating score-variable relationships.
# Length Left Right Bottom Top Diagonal
# [1,] -0.096 0.8785 1.1435 -0.4175 -0.9505 0.5165
# [2,] -0.296 -0.4215 -0.2565 -1.3175 -1.1505 1.2165
# [3,] -0.096 -0.4215 -0.2565 -0.7175 -1.0505 1.7165
# [4,] -0.096 -0.4215 -0.3565 -1.9175 -0.2505 1.5165
# [5,] 0.104 -0.5215 -0.2565 0.9825 -2.9505 1.3165
# [6,] 0.804 0.6785 0.5435 -0.4175 -0.5505 0.9165Project the preprocessed data onto the eigenvectors to obtain scores. Prefer the score matrix from the fitted prcomp() object when available, because it preserves the centering and scaling details.
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] -0.5496481 -0.5063730 -0.27586457 -1.19372812 -1.17050345 0.05836252
# [2,] -2.0186292 -0.6612636 0.50199746 0.01544486 -0.29113718 0.14582450
# [3,] -1.8356755 -1.1753073 -0.04562915 0.18906546 -0.11268124 0.12578824
# [4,] -2.4943672 0.1188916 -0.07652521 0.31613769 -0.08076370 0.07052799
# [5,] -0.7013228 -3.1947667 0.83877387 -0.52104953 0.08262773 0.26879336
# [6,] -0.8458460 -0.4824940 -0.77029691 -1.07530245 -0.09605941 -0.11128017Combine the score matrix and original variables only when they refer to the same rows. Use cor() to calculate associations, and verify that missing observations have not changed pairwise samples unexpectedly.
# 1 2 3 4 5
# 1 1.000000e+00 -7.014633e-18 -4.469949e-16 -7.976627e-16 -1.020254e-15
# 2 -7.014633e-18 1.000000e+00 1.983991e-16 -5.252503e-16 -1.038670e-16
# 3 -4.469949e-16 1.983991e-16 1.000000e+00 4.464784e-16 -1.835115e-16
# 4 -7.976627e-16 -5.252503e-16 4.464784e-16 1.000000e+00 -1.130514e-15
# 5 -1.020254e-15 -1.038670e-16 -1.835115e-16 -1.130514e-15 1.000000e+00
# 6 -6.521526e-16 5.101457e-16 -1.808638e-15 -1.357142e-15 -2.570084e-15
# 6 Length Left Right Bottom Top
# 1 -6.521526e-16 -0.20136050 0.5381321 0.5966697 0.921229423 0.435255426
# 2 5.101457e-16 0.02751049 0.1914237 0.1586673 -0.377020918 0.794217722
# 3 -1.808638e-15 -0.42754733 -0.3538910 -0.4209169 -0.074460283 -0.342054932
# 4 -1.357142e-15 -0.65812392 -0.5566063 -0.4534936 0.056839988 0.247648882
# 5 -2.570084e-15 0.58340637 -0.2804091 -0.3862445 0.020200457 0.037046504
# 6 1.000000e+00 0.04909498 -0.4001151 0.2946144 -0.002898297 -0.008181424
# Diagonal
# 1 -0.870231992
# 2 -0.410109295
# 3 -0.253377217
# 4 0.098959622
# 5 -0.021396513
# 6 -0.007470889Select the columns for the intended variables and components by name where possible. Numeric column positions can be fragile if the data structure changes.
# 1 2
# Length -0.2013605 0.02751049
# Left 0.5381321 0.19142371
# Right 0.5966697 0.15866725
# Bottom 0.9212294 -0.37702092
# Top 0.4352554 0.79421772
# Diagonal -0.8702320 -0.41010930Displaying correlation of original variables with PCs
Set plotting parameters and draw a reference circle only when the selected variable-component coordinates support that interpretation. In a correlation plot, the unit circle marks the maximum absolute correlation.
Label the axes with the components and variance proportions from the fitted PCA. Use the first two component-variable correlations for the plotted positions and keep labels readable without overlap.
par(pty = "s")
ucircle = cbind(cos((0:360)/180*pi), sin((0:360)/180*pi))
plot(ucircle,
type= "l", lty= "solid", col = "blue", lwd = 2,
xlab = "First PC", ylab = "Second PC", main = "Unit variance circle")
abline(h = 0.0, v = 0.0, lty = 2)
text(x = r1,
label = c("Length", "Left", "Right", "Bottom", "Top", "Diagonal"),
cex = 1.2){var m=typeof Object.defineProperties==)
0&&(f=parseInt(a[a.length-1].timestamp,10)+1);var g=null,q=null,G=void 0,H=null,r=(window.location.hash||"#").substring(1),I,J;/^comment-form_/.test(r)?I=r.substring(13):/^c[0-9]+$/.test(r)&&(J=r.substring(1));var V={id:c.postId,data:a,loadNext:function(k){if(f){var h=c.feed+"?alt=json&v=2&orderby=published&reverse=false&max-results=50";f&&(h+="&published-min="+(new Date(f)).toISOString());window.bloggercomments=function(t){f=null;var p=[];if(t&&t.feed&&
t.feed.entry)for(var d,K=0;d=t.feed.entry[K];K++){var l={},e=/blog-(\d+).post-(\d+)/.exec(d.id.$t);l.id=e?e[2]:null;a:{e=void 0;var L=d&&(d.content&&d.content.$t||d.summary&&d.summary.$t)||"";if(d&&d.gd$extendedProperty)for(e in d.gd$extendedProperty)if(d.gd$extendedProperty[e].name=="blogger.contentRemoved"){e=''+L+"";break a}e=L}l.body=e;l.timestamp=Date.parse(d.published.$t)+"";d.author&&d.author.constructor===Array&&(e=d.author[0])&&(l.author={name:e.name?
e.name.$t:void 0,profileUrl:e.uri?e.uri.$t:void 0,avatarUrl:e.gd$image?e.gd$image.src:void 0});d.link&&(d.link[2]&&(l.link=l.permalink=d.link[2].href),d.link[3]&&(e=/.*comments\/default\/(\d+)\?.*/.exec(d.link[3].href))&&e[1]&&(l.parentId=e[1]));l.deleteclass="item-control blog-admin";if(d.gd$extendedProperty)for(var u in d.gd$extendedProperty)d.gd$extendedProperty[u].name=="blogger.itemClass"?l.deleteclass+=" "+d.gd$extendedProperty[u].value:d.gd$extendedProperty[u].name=="blogger.displayTime"&&
(l.displayTime=d.gd$extendedProperty[u].value);p.push(l)}f=p.length<50?null:parseInt(p[p.length-1].timestamp,10)+1;k(p);window.bloggercomments=null};h+="&callback=bloggercomments";var y=P();y.type="text/javascript";O(y,Q(h));document.getElementsByTagName("head")[0].appendChild(y)}},hasMore:function(){return!!f},getMeta:function(k,h){return"iswriter"==k?h.author&&h.author.name==c.authorName&&h.author.profileUrl==c.authorUrl?"true":"":"deletelink"==k?c.baseUri+"/comment/delete/"+c.blogId+"/"+h.id:"deleteclass"==
k?h.deleteclass:""},onReply:function(k,h){g==null&&(g=document.getElementById("comment-editor"),g!=null&&(H=g.style.height,g.style.display="block",q=g.src.split("#")));g&&k&&k!==G&&(document.getElementById(h).insertBefore(g,null),h=q[0]+(k?"&parentID="+k:""),q[1]&&(h=h+"#"+q[1]),g.src=h,g.style.height=H||g.style.height,G=k,g.removeAttribute("data-resized"),g.dispatchEvent(new Event("iframeMoved")))},rendered:!0,initComment:J,initReplyThread:I,config:{maxDepth:c.maxThreadDepth},messages:b};a=function(){if(window.goog&&
window.goog.comments){var k=document.getElementById("comment-holder");window.goog.comments.render(k,V)}};window.goog&&window.goog.comments?a():(window.goog=window.goog||{},window.goog.comments=window.goog.comments||{},window.goog.comments.loadQueue=window.goog.comments.loadQueue||[],window.goog.comments.loadQueue.push(a))},S=["blogger","widgets","blog","initThreadedComments"],T=x,U;S.length&&(U=S.shift());)S.length||R===void 0?T=T[U]&&T[U]!==Object.prototype[U]?T[U]:T[U]={}:T[U]=R;}).call(this);
Thank you for this tutorial. For the final plot, how do I replace the black numbers with values I have in the first column of the data?
ReplyDeleteYou need to convert first variable column to row names of the data set. Then you can add below argument to biplot() function:
Deletepch = rownames(data)
Hope this will work
I truly appreciate the time and work you put into sharing your knowledge. I found this topic to be quite effective and beneficial to me. Thank you very much for sharing. Continue to blog.
ReplyDeleteData Engineering Services
AI & ML Solutions
Data Analytics Services
Data Modernization Services
I truly appreciate the time and work you put into sharing your knowledge. I found this topic to be quite effective and beneficial to me. Thank you very much for sharing. Continue to blog.
ReplyDeleteData Engineering Services
AI & ML Solutions
Data Analytics Services
Data Modernization Services