evals <- c(3.5, 1.0, 0.7, 0.4, 0.25, 0.15)
cumsum(evals) / sum(evals) * 100[1] 58.33333 75.00000 86.66667 93.33333 97.50000 100.00000Activity 3.3 - PCA implementation
SUBMISSION INSTRUCTIONS
- Render to html
- Publish your html to RPubs
- Submit a link to your published solutions
Problem 1
Consider the following 6 eigenvalues from a \(6\times 6\) correlation matrix:
\[\lambda_1 = 3.5, \lambda_2 = 1.0, \lambda_3 = 0.7, \lambda_4 = 0.4, \lambda_5 = 0.25, \lambda_6 = 0.15\]
If you want to retain enough principal components to explain at least 90% of the variability inherent in the data set, how many should you keep?
I should keep 4 because 3 explains only 86.67% variability in the dataset, while 4 explains 93.33% of the data.
Problem 2
The iris data set is a classic data set often used to demonstrate PCA. Each iris in the data set contained a measurement of its sepal length, sepal width, petal length, and petal width. Consider the five irises below, following mean-centering and scaling:
library(tidyverse)
five_irises <- data.frame(
row.names = 1:5,
Sepal.Length = c(0.189, 0.551, -0.415, 0.310, -0.898),
Sepal.Width = c(-1.97, 0.786, 2.62, -0.590, 1.70),
Petal.Length = c(0.137, 1.04, -1.34, 0.534, -1.05),
Petal.Width = c(-0.262, 1.58, -1.31, 0.000875, -1.05)
) %>% as.matrixConsider also the loadings for the first two principal components:
# Create the data frame
pc_loadings <- data.frame(
PC1 = c(0.5210659, -0.2693474, 0.5804131, 0.5648565),
PC2 = c(-0.37741762, -0.92329566, -0.02449161, -0.06694199),
row.names = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
) %>% as.matrixA plot of the first two PC scores for these five irises is shown in the plot below.
Match the ID of each iris (1-5) to the correct letter of its score coordinates on the plot.
(scores <- five_irises %*% pc_loadings) %>% round(3) PC1 PC2
1 0.561 1.762
2 1.572 -1.065
3 -2.440 -2.142
4 0.631 0.415
5 -2.128 -1.1353=a
1=b
4=c
2=d
5=e
Problem 3
These data are taken from the Places Rated Almanac, by Richard Boyer and David Savageau, copyrighted and published by Rand McNally. The nine rating criteria used by Places Rated Almanac are:
- Climate & Terrain
- Housing
- Health Care & Environment
- Crime
- Transportation
- Education
- The Arts
- Recreation
- Economics
For all but two of the above criteria, the higher the score, the better. For Housing and Crime, the lower the score the better. The scores are computed using the following component statistics for each criterion (see the Places Rated Almanac for details):
- Climate & Terrain: very hot and very cold months, seasonal temperature variation, heating- and cooling-degree days, freezing days, zero-degree days, ninety-degree days.
- Housing: utility bills, property taxes, mortgage payments.
- Health Care & Environment: per capita physicians, teaching hospitals, medical schools, cardiac rehabilitation centers, comprehensive cancer treatment centers, hospices, insurance/hospitalization costs index, flouridation of drinking water, air pollution.
- Crime: violent crime rate, property crime rate.
- Transportation: daily commute, public transportation, Interstate highways, air service, passenger rail service.
- Education: pupil/teacher ratio in the public K-12 system, effort index in K-12, accademic options in higher education.
- The Arts: museums, fine arts and public radio stations, public television stations, universities offering a degree or degrees in the arts, symphony orchestras, theatres, opera companies, dance companies, public libraries.
- Recreation: good restaurants, public golf courses, certified lanes for tenpin bowling, movie theatres, zoos, aquariums, family theme parks, sanctioned automobile race tracks, pari-mutuel betting attractions, major- and minor- league professional sports teams, NCAA Division I football and basketball teams, miles of ocean or Great Lakes coastline, inland water, national forests, national parks, or national wildlife refuges, Consolidated Metropolitan Statistical Area access.
- Economics: average household income adjusted for taxes and living costs, income growth, job growth.
In addition to these, latitude and longitude, population and state are also given, but should not be included in the PCA.
Use PCA to identify the major components of variation in the ratings among cities.
places <- read.csv('Data/Places.csv')
head(places) City Climate Housing HlthCare Crime Transp Educ Arts
1 AbileneTX 521 6200 237 923 4031 2757 996
2 AkronOH 575 8138 1656 886 4883 2438 5564
3 AlbanyGA 468 7339 618 970 2531 2560 237
4 Albany-Schenectady-TroyNY 476 7908 1431 610 6883 3399 4655
5 AlbuquerqueNM 659 8393 1853 1483 6558 3026 4496
6 AlexandriaLA 520 5819 640 727 2444 2972 334
Recreat Econ Long Lat Pop
1 1405 7633 -99.6890 32.5590 110932
2 2632 4350 -81.5180 41.0850 660328
3 859 5250 -84.1580 31.5750 112402
4 1617 5864 -73.7983 42.7327 835880
5 2612 5727 -106.6500 35.0830 419700
6 1018 5254 -92.4530 31.3020 135282A.
Perform PCA on the data set, using columnsClimate through Econ, and create a scree plot. Interpret your findings from this plot. What percent of total variability is explained by the first two PCs?
library(FactoMineR)
library(factoextra)Welcome to factoextra!Want to learn more? See two factoextra-related books at https://www.datanovia.com/library/principal-component-methodsplaces_pca <- places %>%
select(Climate:Econ) %>%
PCA(ncp = 9, graph = FALSE)
places_pca$eig %>% round(3) eigenvalue percentage of variance cumulative percentage of variance
comp 1 3.408 37.870 37.870
comp 2 1.214 13.489 51.359
comp 3 1.141 12.683 64.042
comp 4 0.921 10.232 74.274
comp 5 0.753 8.370 82.644
comp 6 0.631 7.006 89.650
comp 7 0.493 5.478 95.128
comp 8 0.318 3.534 98.662
comp 9 0.120 1.338 100.000The first two PCs explain 37.9% + 13.5% = 51.4% of the total variability, a 2D picture captures only about half of the information in the nine ratings.
B.
Interpret the first two principal components by examining the loadings (plot(s) of the loadings may be helpful). Which variables will be used to separate cities along the first and second principal axes, and how? Make sure to discuss the signs of the loadings, not just their contributions!
library(patchwork)
c1 <- fviz_contrib(places_pca, choice = 'var', axes = 1) +
scale_x_discrete(limits = rev, name = '') +
coord_flip() +
theme_classic(base_size = 14) +
labs(title = 'Contribution to PC1')
c2 <- fviz_contrib(places_pca, choice = 'var', axes = 2) +
scale_x_discrete(limits = rev, name = '') +
coord_flip() +
theme_classic(base_size = 14) +
labs(title = 'Contribution to PC2')
c1 + c2
PC1: Arts and Health Care (about 21% each) dominate, and Housing and Transportation also cross the line. Recreation is just under it. Econ barely matters.
PC2: Education and Econ (about 22–23%) dominate, followed by Recreation and Crime.
# Signed loadings
loadings_df <- data.frame(places_pca$var$coord) %>%
rownames_to_column('Variable')
base_plot <- ggplot(loadings_df) +
theme_classic(base_size = 14) +
labs(y = '', x = 'Loading') +
geom_vline(aes(xintercept = 0), linetype = 2)
l1 <- base_plot +
geom_col(aes(y = reorder(Variable, Dim.1), x = Dim.1)) +
labs(title = 'Dimension 1')
l2 <- base_plot +
geom_col(aes(y = reorder(Variable, Dim.2), x = Dim.2)) +
labs(title = 'Dimension 2')
l1 + l2
Dimension 1: Every bar is positive. A city with a high PC1 score is high on all nine ratings, mostly Arts, Health Care, Housing, and Transportation. Since high Housing and Crime are bad, these cities are also expensive and have more crime, which tend to be bigger cities.
Dimension 2: Econ, Recreation, Crime, Housing, and Climate point positive. Education, Health Care, Arts, and Transportation are negative. A city with a high PC2 score has a strong economy and recreation, but weaker education and health care. A city with a low PC2 score is the opposite.
# Loadings plot (arrows)
fviz_pca_var(places_pca, axes = c(1, 2), repel = TRUE)
All arrows point to the right, so every variable is positively correlated with PC1. Arts, Health Care, and Transportation point right and slightly down. Recreation, Crime, and Housing cluster together, showing they’re positively correlated with each other, while the short Climate arrow shows it isn’t captured well by the first two PCs.
C.
Add the first two PC scores to the places data set. Create a biplot of the first 2 PCs, using repelled labeling to identify the cities. Which are the outlying cities and what characteristics make them unique?
library(ggrepel)
# Add first two PC scores to the data set
places_update <- bind_cols(places, places_pca$ind$coord[, 1:2])
# Most extreme cities
places_update %>%
arrange(desc(Dim.1)) %>%
select(City, Dim.1, Dim.2) %>%
head(10) City Dim.1 Dim.2
1 New-YorkNY 12.445179 -2.0648866
2 San-FranciscoCA 7.402662 1.2288687
3 Los-AngelesLong-BeachCA 7.257344 0.4778784
4 ChicagoIL 6.474759 -3.0918369
5 BostonMA 6.310655 -1.6111921
6 WashingtonDC-MD-VA 6.195895 -2.2683238
7 PhiladelphiaPA-NJ 4.772628 -3.0779683
8 BaltimoreMD 4.303905 -0.9694944
9 NewarkNJ 4.139065 -0.7909262
10 SeattleWA 3.990213 1.5896876#| fig-alt: "Biplot of first two PCs for Places Rated, with City Labels"
#| fig-width: 10
#| fig-height: 8
#| warning: false
fviz_pca(places_pca, axes = c(1, 2),
label = 'var', # only label the arrows; cities labeled below
col.ind = 'grey50',
repel = TRUE) +
geom_text_repel(aes(x = Dim.1, y = Dim.2, label = City),
data = places_update,
size = 2.5,
max.overlaps = 10,
show.legend = FALSE) +
theme_classic(base_size = 14)
New York is the biggest outliar, far to the right on PC1 because it has the highest Arts, Health Care, Transportation, and Crime ratings and very expensive housing. San Francisco, Los Angeles, Boston, Washington DC, Chicago, and Philadelphia also sit far right as large,rich, costly cities. San Francisco and LA are above 0 on PC2 because of their climate and recreation, while Chicago and Philadelphia sit below 0 because of strong education and health care relative to economic growth. At the top of the plot, Las Vegas, Atlantic City, and Midland TX stand out for strong economies,recreation, weaker education, and health care.
Problem 4
The data we will look at here come from a study of malignant and benign breast cancer cells using fine needle aspiration conducted at the University of Wisconsin-Madison. The goal was determine if malignancy of a tumor could be established by using shape characteristics of cells obtained via fine needle aspiration (FNA) and digitized scanning of the cells.
The variables in the data file you will be using are:
- ID - patient identification number (not used in PCA)
- Diagnosis determined by biopsy - B = benign or M = malignant
- Radius: mean of distances from center to points on the perimeter
- Texture: standard deviation of gray-scale values
- Smoothness: local variation in radius lengths
- Compactness: perimeter^2 / area - 1.0
- Concavity: severity of concave portions of the contour
- Concavepts: number of concave portions of the contour
- Symmetry: measure of symmetry of the cell nucleus
- FracDim: fractal dimension; “coastline approximation” - 1
bc_cells <- read.csv('Data/BreastDiag.csv')
head(bc_cells) Diagnosis Radius Texture Smoothness Compactness Concavity ConcavePts Symmetry
1 M 17.99 10.38 0.11840 0.27760 0.3001 0.14710 0.2419
2 M 20.57 17.77 0.08474 0.07864 0.0869 0.07017 0.1812
3 M 19.69 21.25 0.10960 0.15990 0.1974 0.12790 0.2069
4 M 11.42 20.38 0.14250 0.28390 0.2414 0.10520 0.2597
5 M 20.29 14.34 0.10030 0.13280 0.1980 0.10430 0.1809
6 M 12.45 15.70 0.12780 0.17000 0.1578 0.08089 0.2087
FracDim
1 0.07871
2 0.05667
3 0.05999
4 0.09744
5 0.05883
6 0.07613A.
My analysis suggests 3 PCs should be retained. Support or refute this suggestion. What percent of variability is explained by the first 3 PCs?
bc_pca <- PCA(bc_cells, ncp = 8, quali.sup = 'Diagnosis', graph = FALSE)
bc_pca$eig %>% round(3) eigenvalue percentage of variance cumulative percentage of variance
comp 1 4.287 53.589 53.589
comp 2 1.823 22.794 76.383
comp 3 0.826 10.322 86.704
comp 4 0.499 6.233 92.937
comp 5 0.372 4.654 97.591
comp 6 0.092 1.152 98.743
comp 7 0.069 0.860 99.602
comp 8 0.032 0.398 100.000#| fig-alt: "Scree plot for breast cancer cell PCA"
fviz_eig(bc_pca, geom = 'line') +
theme_classic(base_size = 16)
The first 3 PCs explain 86.7% of the total variability (53.6% + 22.8% + 10.3%), and the scree plot levels off after PC3, with PC4 adding only 6.2%. 3 makes sense.
B.
Interpret the first 3 principal components by examining the eigenvectors/loadings. Discuss.
contrib_plot <- function(ax) {
fviz_contrib(bc_pca, choice = 'var', axes = ax) +
scale_x_discrete(limits = rev, name = '') +
coord_flip() +
theme_classic(base_size = 12) +
labs(title = paste0('Contribution to PC', ax))
}
contrib_plot(1) + contrib_plot(2) + contrib_plot(3)
#| fig-width: 12
#| fig-height: 4.5
bc_loadings <- data.frame(bc_pca$var$coord) %>%
rownames_to_column('Variable')
loading_plot <- function(dim_col, ttl) {
ggplot(bc_loadings, aes(y = reorder(Variable, .data[[dim_col]]),
x = .data[[dim_col]])) +
geom_col(fill = 'steelblue') +
geom_vline(xintercept = 0, linetype = 2) +
scale_x_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
theme_classic(base_size = 12) +
labs(y = '', x = 'Loading', title = ttl)
}
loading_plot('Dim.1', 'PC1') + loading_plot('Dim.2', 'PC2') + loading_plot('Dim.3', 'PC3')
bc_pca$eig[1:3, ] %>% round(3) eigenvalue percentage of variance cumulative percentage of variance
comp 1 4.287 53.589 53.589
comp 2 1.823 22.794 76.383
comp 3 0.826 10.322 86.704PC1 (53.6%) is an overall shape-irregularity score. Every loading is positive, led by Compactness, Concavity, and Concave Points, so high scores mean large irregular cells with many deep indentations. PC2 (22.8%) contrasts boundary complexity with size. FracDim, Smoothness, and Symmetry load positively while Radius and Texture load negatively, so high scores mean small cells with jagged borders. PC3 (10.3%) is almost entirely Texture (+0.82), separating cells by how much their gray-scale values are different.
C.
Examine a biplot of the first two PCs. Incorporate the third PC by sizing the points by this variable. (Hint: you might want to add all 3 PC scores to the data set. Use fviz_pca to set up a biplot, but set col.ind='white'. Then use geom_point() to maintain full control over the point mapping.) Color-code by whether the cells are benign or malignant.
bc_update <- bind_cols(bc_cells, bc_pca$ind$coord[, 1:3])
#| fig-alt: "Biplot of PC1 vs PC2, colored by diagnosis and sized by PC3"
#| fig-width: 9
#| fig-height: 7
fviz_pca(bc_pca, axes = c(1, 2),
col.ind = 'white',
label = 'var',
repel = TRUE) +
geom_point(aes(x = Dim.1, y = Dim.2, color = Diagnosis, size = Dim.3),
data = bc_update, alpha = 0.5) +
scale_color_manual(values = c(B = 'steelblue', M = 'firebrick'),
labels = c('Benign', 'Malignant')) +
scale_size_continuous(range = c(0.5, 4)) +
labs(color = 'Diagnosis', size = 'PC3 score', title = '') +
theme_classic(base_size = 14)
bc_update %>%
group_by(Diagnosis) %>%
summarize(across(Dim.1:Dim.3, mean))# A tibble: 2 × 4
Diagnosis Dim.1 Dim.2 Dim.3
<chr> <dbl> <dbl> <dbl>
1 B -1.11 0.426 -0.00172
2 M 1.87 -0.717 0.00289#| fig-width: 10
#| fig-height: 4
bc_update %>%
pivot_longer(Dim.1:Dim.3, names_to = 'PC', values_to = 'Score') %>%
ggplot(aes(x = Diagnosis, y = Score, fill = Diagnosis)) +
geom_boxplot() +
facet_wrap(~PC) +
scale_fill_manual(values = c(B = 'steelblue', M = 'firebrick')) +
theme_classic(base_size = 14)
Answer the following:
What characteristics distinguish malignant from benign cells?
Malignant cells (red) fall mostly on the right side of the bi plot, in the direction of the Compactness, Concavity, Concave Points, and Radius arrows, while benign cells (blue) cluster on the left. So malignant cells tend to be larger, with more irregular, concave, and compact borders, and benign cells are smaller, smoother, and rounder. Malignant cells also sit lower on PC2, toward the Radius and Texture arrows, which again reflects their larger size and more uneven texture.
Of the 3 PCs, which does the best job of differentiating malignant from benign cells?
PC1. The group means are far apart and the box plots barely overlap. PC2 shows a smaller shift with considerable overlap. PC3 does pretty much nothing, both groups average about 0, the box plots are nearly identical, and large and small points are mixed throughout both colors in the bi plot.