Activity 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?

PC1: 3.5 / 6 = 58.3% PC2: (3.5 + 1.0) / 6 = 75.0% PC3: (3.5 + 1.0 + 0.7) / 6 = 86.7% PC4: (3.5 + 1.0 + 0.7 + 0.4) / 6 = 93.3%

To reach at least 90%, you need to keep 4 principal components.

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.matrix

Consider 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.matrix

A 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.

Iris 1 → b Iris 2 → d Iris 3 → a Iris 4 → c Iris 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 135282

A.

If you want to explore this data set in lower dimensional space using the first \(k\) principal components, how many would you use, and what percent of the total variability would these retained PCs explain? Use a scree plot to help you answer this question.

places <- read.csv("Data/Places.csv")
rating_cols <- c("Climate", "Housing", "HlthCare", "Crime",
                 "Transp", "Educ", "Arts", "Recreat", "Econ")

X <- places[, rating_cols]
X_scaled <- scale(X)
pca_out <- prcomp(X_scaled, center = TRUE, scale. = TRUE)

plot(pca_out, type = "l", main = "Scree Plot: Places Rated PCA")

var_exp <- pca_out$sdev^2 / sum(pca_out$sdev^2)
var_exp

[1] 0.37869909 0.13488624 0.12683102 0.10232420 0.08369832 0.07006243 0.05478308
[8] 0.03533761 0.01337801

cum_var <- cumsum(var_exp)
cum_var

[1] 0.3786991 0.5135853 0.6404163 0.7427405 0.8264389 0.8965013 0.9512844
[8] 0.9866220 1.0000000

PC1 ≈ 37.9% PC2 ≈ 13.5% PC3 ≈ 12.7% So the first three PCs explain roughly 64% of the total variability, which is a reasonable amount for exploring the cities in a lower-dimensional (2D–3D) space.

B.

Interpret the retained 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!

# Loadings (PC1 and PC2)
loadings <- pca_out$rotation[, 1:2]
loadings

               PC1        PC2
Climate  0.2064140  0.2178353
Housing  0.3565216  0.2506240
HlthCare 0.4602146 -0.2994653
Crime    0.2812984  0.3553423
Transp   0.3511508 -0.1796045
Educ     0.2752926 -0.4833821
Arts     0.4630545 -0.1947899
Recreat  0.3278879  0.3844746
Econ     0.1354123  0.4712833

biplot(pca_out, choices = 1:2, cex = 0.7, main = "PCA Loadings Plot")

PC1 All loadings positive. Represents overall amenities and services. High PC1 = stronger arts, health care, transportation, recreation, and housing. Low PC1 = weaker across those same areas.

PC2 Mix of positive and negative loadings. Separates cities based on a contrast between: Economic/Recreation strengths (positive) Education/HealthCare/Arts strengths (negative) High vs low PC2 reflects different profiles of city strengths rather than overall quality.

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(ggplot2)
library(ggrepel)

places$PC1 <- pca_out$x[,1]
places$PC2 <- pca_out$x[,2]

head(places[, c("City", "PC1", "PC2")])

                       City        PC1         PC2
1                 AbileneTX -1.0401799  0.89376897
2                   AkronOH  0.4398136  0.07506618
3                  AlbanyGA -1.8755393  0.06979169
4 Albany-Schenectady-TroyNY  0.9107414 -1.81758215
5             AlbuquerqueNM  2.1492475  0.32885808
6              AlexandriaLA -1.7879611 -0.78120167

load_df <- data.frame(
  variable = rownames(pca_out$rotation),
  PC1 = pca_out$rotation[,1],
  PC2 = pca_out$rotation[,2]
)

load_df$PC1s <- load_df$PC1 * 3
load_df$PC2s <- load_df$PC2 * 3

# Biplot 
ggplot() +
  geom_point(data = places, aes(PC1, PC2), color="steelblue", size=2) +
  geom_text_repel(data = places, aes(PC1, PC2, label = City), size = 3) +
  
  # arrows
  geom_segment(data = load_df,
               aes(x = 0, y = 0, xend = PC1s, yend = PC2s),
               arrow = arrow(length = unit(0.25, "cm")),
               color = "firebrick") +
  geom_text_repel(data = load_df,
                  aes(PC1s, PC2s, label = variable),
                  color = "firebrick", size = 3.2) +

  theme_minimal() +
  labs(title = "Biplot of the First Two Principal Components")

Warning: ggrepel: 304 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

Cities that often stand out strongly on PC1 or PC2 include: New York, San Francisco, Los Angeles

Why they’re outliers: Extremely strong scores in Arts, Economics, Transportation High population density and very unique housing/commute characteristics They tend to lie far right on PC1 or high on PC2.

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.07613

A.

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_cells <- read.csv("Data/BreastDiag.csv")
vars <- c("Radius","Texture","Smoothness","Compactness",
          "Concavity","ConcavePts","Symmetry","FracDim")

X <- bc_cells[, vars]
X_scaled <- scale(X)
pca_bc <- prcomp(X_scaled, center = TRUE, scale. = TRUE)

# Scree plot
plot(pca_bc, type = "l", main = "Scree Plot: Breast Cancer PCA")

var_exp <- pca_bc$sdev^2 / sum(pca_bc$sdev^2)
var_exp

[1] 0.535890830 0.227935572 0.103215568 0.062328991 0.046536586 0.011517961
[7] 0.008597524 0.003976969

cum_var <- cumsum(var_exp)
cum_var

[1] 0.5358908 0.7638264 0.8670420 0.9293710 0.9759075 0.9874255 0.9960230
[8] 1.0000000

The scree supports the decision to retain 3 PCs. The first three principal components together explain more than half of the total variability in the breast cancer cell measurements, which is sufficient to summarize most of the information in the dataset.

B.

Interpret the first 3 principal components by examining the eigenvectors/loadings. Discuss.

PC1 represents overall size and irregularity of the nucleus. It loads heavily on radius, texture, compactness, concavity, and concave points. High PC1 indicates large, highly irregular cells typical of malignant tumors.

PC2 contrasts shape symmetry vs. texture. High PC2 corresponds to smoother, more symmetrical nuclei, while low PC2 reflects more textural variation and less symmetry.

PC3 captures a finer distinction between smooth vs. deeply concave nuclear borders. High PC3 indicates smooth edges; low PC3 indicates strong concavity and jagged boundaries.

C.

Examine a biplot of the first two PCs. Incorporate the third PC by sizing the points by this variable. (Hint: 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. Answer the following:

What characteristics distinguish malignant from benign cells?
Of the 3 PCs, which does the best job of differentiating malignant from benign cells?

library(ggplot2)
library(ggrepel)

# Scores (PC1, PC2, PC3)
scores <- data.frame(
  PC1 = pca_bc$x[,1],
  PC2 = pca_bc$x[,2],
  PC3 = pca_bc$x[,3],
  Diagnosis = bc_cells$Diagnosis
)

loadings <- data.frame(
  Variable = rownames(pca_bc$rotation),
  PC1 = pca_bc$rotation[,1],
  PC2 = pca_bc$rotation[,2]
)

loadings$PC1s <- loadings$PC1 * 4
loadings$PC2s <- loadings$PC2 * 4

# Biplot
ggplot() +
  geom_point(data = scores,
             aes(x = PC1, y = PC2, color = Diagnosis, size = PC3),
             alpha = 0.7) +
  geom_segment(data = loadings,
               aes(x = 0, y = 0, xend = PC1s, yend = PC2s),
               arrow = arrow(length = unit(0.25, "cm")),
               color = "black") +
  geom_text_repel(data = loadings,
                  aes(x = PC1s, y = PC2s, label = Variable),
                  size = 3.2) +
  
  scale_color_manual(values = c("B" = "steelblue", "M" = "red")) +
  labs(
    title = "Breast Cancer PCA Biplot (PC1 vs PC2, Point Size = PC3)",
    x = "PC1",
    y = "PC2"
  ) +
  theme_minimal()

Malignant = large + irregular + concave Benign = small + smooth + symmetric.

PC1 is the primary component that separates malignant from benign cells, because it captures differences in overall size, concavity, compactness, and irregularity of the cell nuclei.