MouthWidth TotalWidth TotalHeight BaseWidth StemWidth StemHeight
1 13 21 23 14 7 8
2 14 14 24 19 5 9
3 19 23 24 20 6 12
4 17 18 16 16 11 8
5 19 20 16 16 10 7
6 12 20 24 17 6 93.3 - PCA implementation
6D example: Prehistoric goblets
Suppose you have been asked to analyze data on 25 prehistoric goblets from Thailand (Professor C.F.W. Higham, University of Otago, as taken from Manly, B.F.J. 1986. Multivariate Statistical Methods: A Primer. Chapman and Hall, London, 159pp.)
Let’s consider three of these:
- \(X_1\) = mouth width
- \(X_2\) = total width
- \(X_3\) = total height
- \(X_4\) = base width
- \(X_5\) = stem width
- \(X_6\) = stem height
Reading and plotting
- Recall questions:
- These variables are correlated. Do we need all 6 to characterize the goblets?
- Is there a linear combination of these variables that reduces these variables to 2 dimensions that works well to characterize these goblets?
Scree plot with fviz_eig()
- A scree plot is a common plot in PCA that visualizing the variability explained by each PCA.
- Scree is a geological term for a collection of broken rock fragments at the base of steep rocky mass:
Image source: Wikipedia
Scree plot interpretation
- Look for the spot where the cliff ends and the scree begins (point of diminishing returns)
- Also called “elbow” or “kink” in the plot
- These are the number of PCs that adequately approximate the full-dimensional data
- Any PCs beyond the sharp elbow add little additional information
Interpretation:
- We see the cliff ends at PC 2
- Scree begins at PC 3
- 2 dimensions is sufficient for representing this 6D data set
Scree plot options
- Scree plots are often visualized just with a line for a simpler look
- Sometimes represent the eigenvalues instead of the % variability explained
- Both versions easy to get with
fviz_eig(); note also the seamlessggplotintegration:
Variable contributions with fviz_contrib()
- “Variable contributions” to each PC are measured by the squared loadings
fviz_contrib()plots the squared loadings, which must sum to 1 (remember, each vector of loadings is normalized)- A reference line is placed at \(100\times 1/p\), 16.67% for us: the contribution of each variable if each contributed equally to that PC.
- Variables above this threshold are weighted more heavily in that direction than variables below.
Interpreting variable contribution plot: 1st PC
fviz_contrib(goblets_pca, choice = 'var', axes = 1) +
theme_classic(base_size = 16) +
labs(x = 'Variable',
title = 'Contribution to first dimension')
BaseWidth TotalWidth StemHeight TotalHeight MouthWidth StemWidth
0.4618586 0.4515367 0.4381125 0.4111609 0.3660233 0.2963653 BaseWidth TotalWidth StemHeight TotalHeight MouthWidth StemWidth
0.21331336 0.20388536 0.19194254 0.16905330 0.13397304 0.08783241 [1] 1- Base width, total width, stem height, and total height contribute the most to the first PC.
- Mouth width and stem width are less for separating goblets in the first PC direction
- As all loadings are \(>0\), generally large goblets will have large PC1 values and small goblets will have small PC1 values
Interpreting variable contribution plot: 2nd PC
fviz_contrib(goblets_pca, choice = 'var', axes = 2) +
theme_classic(base_size = 16)+
labs(x = 'Variable',
title = 'Contribution to second dimension')
StemWidth MouthWidth TotalWidth BaseWidth StemHeight TotalHeight
0.68277080 0.48592912 -0.03412653 -0.11457532 -0.29768029 -0.44135161 StemWidth MouthWidth TotalHeight StemHeight BaseWidth TotalWidth
0.46617596 0.23612711 0.19479125 0.08861356 0.01312750 0.00116462 [1] 1- Stem width and mouth width (least important for the first PC) are now the two most important variables for separating goblets in the second PC direction.
- Total height is third most important, but has a negative loading
- Goblets with large stem and mouth widths, and small heights (i.e., short, squat goblets) will have large PC2 values
Loadings plot with fviz_pca_var()
- While the variable contribution plot shows us importance, it doesn’t convey how the variables explain variability in the first 2 PC directions.
- That’s the job of the loadings plot, which plots the loadings as arrow vectors.
- Total width and base width point almost directly horizontal, as we would expect given the variable contribution plots: the first PC is determined primarily by these variables
- Stem width and mouth width point more in the vertical direction than any other variable, consistent with these being the highest contributors to PC2.
- All arrows point in somewhat similar directions, consistent with notable positive correlation between all of these variables
Biplot with fviz_pca()
- We can put this all together, adding the arrows on top of a plot of the first two PC scores with
fviz_pca() - One point = one goblet’s first two PC scores
- One arrow = one variable’s loading
PC1 PC2
[1,] 0.4724326 -0.025587670
[2,] 0.3551757 -0.747416633
[3,] 2.5682176 -0.005956072
[4,] 0.9455738 2.363098749
[5,] 1.1203843 2.448935880
[6,] 0.6508633 -0.746100727
[7,] 0.4394648 -0.660112674
[8,] 0.5621160 -0.273163556
[9,] -1.6920335 0.185212389
[10,] -2.1054544 0.830051502
[11,] 1.1112325 -1.447189946
[12,] 0.8585801 0.579560252
[13,] -1.3213315 -0.237320609
[14,] 1.3930101 -0.523917666
[15,] 1.1522576 -0.212222372
[16,] 0.5217002 -0.535724932
[17,] -0.0994333 1.617947532
[18,] 1.7763959 1.549731588
[19,] 0.9515815 -0.642113521
[20,] 2.2715217 -0.657717663
[21,] 0.9755721 -0.802412037
[22,] -4.0334505 0.085097503
[23,] -5.0994230 -0.333415717
[24,] -4.9173947 -0.224676863
[25,] 1.1424410 -1.584586738Interpreting biplot
- Goblets like 3 and 20 have very wide Total and Base widths
- Goblets like 23 and 24 have very narrow Total and Base widths
- Goblet 20 is much taller than Goblet 10
- Goblets like 4 and 5 have similar Total and Base width to goblets like 15 and 14, but have much widers Stems and Mouths.
Pairs plot
- What kind of correlations are we dealing with?
- calories, fat (both total and saturated) and protein all correlated with each other
- Cholesteral highly correlated with protein and fat
- Sugar, fiber, carbs, and cholesterol are kind of “on their own”
How many PCs?
Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
Standard deviation 2.2601 1.1632 0.93582 0.91378 0.61020 0.50978 0.34841 0.26973 0.04262
Proportion of Variance 0.5676 0.1503 0.09731 0.09278 0.04137 0.02888 0.01349 0.00808 0.00020
Cumulative Proportion 0.5676 0.7179 0.81520 0.90798 0.94935 0.97823 0.99171 0.99980 1.00000
- Could argue the “scree” starts at 5 and above
- Consider retaining first 4 PCs to approximate data? (Explains 90.8% of total variability)
What’s contributing to the PCs?
PC1 PC2 PC3 PC4
Calories 0.42671188 -0.06019828 -0.04565786 -0.23226760
TotalFat 0.41555181 0.10272038 0.02777273 -0.19041622
SaturatedFat 0.39977015 0.07660561 0.22907804 0.08035464
Cholesterol 0.35766578 0.19339534 0.16139158 0.32754323
Sodium 0.36317299 0.10118529 -0.23634847 -0.10139831
TotalCarbs 0.18874417 -0.59221814 -0.24040894 -0.54335240
Fiber 0.13019795 -0.46002736 -0.53815934 0.65806472
Sugar 0.06708248 -0.59171566 0.71730148 0.16443492
Protein 0.41051448 0.14085855 0.02874194 0.17450080- Note the “equal contribution” threshold = \(100\times 1/10 =10\%\)
- Calories, fat, and protein are driving PC1
- PC2 is driven by carbs, sugar and fiber. Sugar and fiber are main drivers of PC3.
- Cholesterol third highest contributor to PC4.
Sign matters!
- It takes a bit more work to plot the raw loadings (must create out own data frame), but in this way we can see their signs which provides insight the squared loadings don’t.
loadings_df <- data.frame(food_pca$rotation) %>%
rownames_to_column('Variable')
base_plot <- ggplot(loadings_df) +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45,hjust=1,vjust=1)) +
labs(x = 'Variable') +
geom_hline(aes(yintercept= 0), linetype=2)
l1 <- base_plot +
geom_col(aes(x = reorder(Variable,PC1, decreasing=TRUE), y = PC1))
l2 <- base_plot+
geom_col(aes(x = reorder(Variable,PC2, decreasing=TRUE), y = PC2))
l3 <- base_plot +
geom_col(aes(x = reorder(Variable,PC3, decreasing=TRUE), y = PC3))
l4 <- base_plot +
geom_col(aes(x = reorder(Variable,PC4, decreasing=TRUE), y = PC4))
(l1+l2)/(l3+l4) +
plot_annotation(title = 'Loadings plots for fast food PCA',
theme=theme(plot.title = element_text(hjust = 0.5)))- Foods low in fiber, sugar, carbs will have high PC2 scores.
- PC3 contrasts foods high in sugar and low in fiber.
- PC4 contrasts foods high in fiber and low in carbs.
PCA biplot + bells and whistles
- We can easily incorporate more information (ggrepeled-point labeling, color-coding) in PCA biplots
- Default graph: ugh!!
Cleaning it up
- Here we are also suppressing the point labeling, labeling only the variable arrows (even repelling doesn’t save us if we label everything!), and cleaning up the legend
Interpretation: axes 1 and 2
- The “big points” are the mean vectors of the PCs for each food type
- Protein, sodium, saturated fat, total fat, cholesterol all highly correlated
- Food items that are high in these are obviously sandwiches
- Desserts are low in these
- These are not that correlated with carbs, sugar, fiber
- Low in sugar/fiber/carbs include chicken and wraps
Plot result
- Notice how PCs 3 and 4 are much more focused on the sugar, cholesterol, fiber, and carb content - the “low-correlation” variables in our initial pairs plot
- PCs 3 and 4 highlight variability in the salads, desserts, and sides - notice how clumped in the middle the sandwiches and chicken are in the PC3v4 plot, when they were the big story in the PC1v2 plot
- Notice how PCs 3 and 4 contrast fiber, sugar, and carbs, which we knew from the loadings plots)