Activity 3.3 - PCA implementation

SUBMISSION INSTRUCTIONS

1. Render to html
2. Publish your html to RPubs
3. 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?

Calculation	Proportion	Cumulative
3.5 / 6	.583	58.3%
1 / 6	.167	75.0%
0.7 / 6	.117	86.7%
0.4 / 6	.067	93.3%

Having 4 principal components explains 93.3% of variability, so 4 is the minimum number we should keep.

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.

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(tidyverse)
library(FactoMineR)
library(factoextra)

Welcome to factoextra!

Want to learn more? See two factoextra-related books at https://www.datanovia.com/library/principal-component-methods

places_pca <- (places %>%
                 select(Climate:Econ) %>% 
                 PCA(ncp = 9, graph = FALSE))

places_pca$eig

       eigenvalue percentage of variance cumulative percentage of variance
comp 1  3.4082918              37.869909                          37.86991
comp 2  1.2139762              13.488624                          51.35853
comp 3  1.1414791              12.683102                          64.04163
comp 4  0.9209178              10.232420                          74.27405
comp 5  0.7532849               8.369832                          82.64389
comp 6  0.6305619               7.006243                          89.65013
comp 7  0.4930477               5.478308                          95.12844
comp 8  0.3180385               3.533761                          98.66220
comp 9  0.1204021               1.337801                         100.00000

(fviz_eig(places_pca) +
    labs(title = 'Scree plot') +
    theme_classic(base_size = 14))

Cliff ends at PC1, scree begins at PC2 therefore, 2 dimensions is not sufficient for representing this 9D data set.

The first two PCs explain 37.87% + 13.49% = 51.36% of total variability in 9-dimensional space

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)

(places_pca$var$coord[, 1:2] %>%
    round(3))

         Dim.1  Dim.2
Climate  0.381  0.240
Housing  0.658  0.276
HlthCare 0.850 -0.330
Crime    0.519  0.392
Transp   0.648 -0.198
Educ     0.508 -0.533
Arts     0.855 -0.215
Recreat  0.605  0.424
Econ     0.250  0.519

c1 <- (fviz_contrib(places_pca, choice = 'var', axes = 1) +
         theme_classic(base_size = 14) +
         labs(x = 'Variable', title = 'cont. to first dimension'))
c2 <- (fviz_contrib(places_pca, choice = 'var', axes = 2) +
         theme_classic(base_size = 14) +
         labs(x = 'Variable', title = 'cont. to second dimension'))

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'))

(c1 + c2) / (l1 + l2)

fviz_pca_var(places_pca, axes = c(1, 2), repel = TRUE)

PC1: All loadings are positive, so PC1 is an overall size axis. Arts and HlthCare contribute the most, followed by Housing, Transp, and Recreat; Econ and Climate contribute little. We can intrepret this as cities with large PC1 scores have arts, health care, transportation, and recreation, but also expensive housing and more crime. Cities with small PC1 scores are smaller places with fewer amenities but cheaper housing and less crime.

PC2: PC2 is driven mainly by Econ (positive) and Educ (negative). Recreat, Crime, Housing, and Climate also load positively, while HlthCare, Arts, and Transp load negatively. Cities with large PC2 scores have strong economies and recreation but higher crime and housing costs. Cities with small PC2 scores are strong in education and health care relative to their economies, like university or medical-center towns.

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)

places_update <- bind_cols(places, places_pca$ind$coord[, 1:2])

(fviz_pca(places_pca, axes = c(1, 2),
          label = 'var',
          repel = TRUE,
          col.ind = 'gray50',
          pointsize = 1.5) +
    geom_text_repel(aes(x = Dim.1, y = Dim.2, label = City),
                    data = places_update,
                    size = 3,
                    max.overlaps = 5) +
    labs(title = 'Biplot of first two PCs'))

(places_update %>%
    select(City, Dim.1, Dim.2) %>%
    arrange(desc(abs(Dim.1))) %>%
    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

(places_update %>%
    select(City, Dim.1, Dim.2) %>%
    arrange(desc(abs(Dim.2))) %>%
    head(10))

                                        City      Dim.1     Dim.2
1                                Las-VegasNV  1.3559807  3.344152
2                                  ChicagoIL  6.4747591 -3.091837
3                          PhiladelphiaPA-NJ  4.7726283 -3.077968
4                            Atlantic-CityNJ  1.3712239  3.053620
5                                  MidlandTX -0.1099815  3.034607
6                 Salinas-Seaside-MontereyCA  1.0832172  2.877995
7                                LafayetteLA -0.2355063  2.545723
8                            Miami-HialeahFL  3.6513637  2.545412
9  West-Palm-Beach-Boca-Raton-Delray-BeachFL  2.2059328  2.498020
10                          Oxnard-VenturaCA  0.3094478  2.472838

Some of the nations biggest cities are outliers, specifically, New York appears to be an outlier city with a PC1 score around 12, almost twice that of the next-closest cities. That means it rates far above average on most of the heavily weighted PC1 variables at once: Arts, HlthCare, Transp, and Recreat. Other big cities like San Francisco and Los Angeles share these traits, but not as extreme as New York.

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_pca <- (bc_cells %>%
             select(Radius:FracDim) %>%
             PCA(ncp = 8, graph = FALSE))

bc_pca$eig

       eigenvalue percentage of variance cumulative percentage of variance
comp 1 4.28712664             53.5890830                          53.58908
comp 2 1.82348458             22.7935572                          76.38264
comp 3 0.82572454             10.3215568                          86.70420
comp 4 0.49863192              6.2328991                          92.93710
comp 5 0.37229268              4.6536586                          97.59075
comp 6 0.09214368              1.1517961                          98.74255
comp 7 0.06878019              0.8597524                          99.60230
comp 8 0.03181576              0.3976969                         100.00000

(fviz_eig(bc_pca) +
    labs(title = 'Scree plot') +
    theme_classic(base_size = 14))

bc_pca$eig[3, 'cumulative percentage of variance']

[1] 86.7042

I would keep 3 PCs, which together explain 86.70% of the total variability in the 8 shape variables.

The scree plot shows a pretty clear drop after the first three PCs. PC1 explains 53.6%, PC2 explains 22.8%, and PC3 explains 10.3%. After that, the plot starts to flatten out, with PC4 only adding 6.2% and PCs 6 through 8 each adding less than 1.2%. Because the remaining PCs don’t add much information, I think keeping 3 PCs is enough to represent the original 8 variables without losing too much information.

B.

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

(bc_pca$var$coord[, 1:3] %>%
    round(3))

            Dim.1  Dim.2  Dim.3
Radius      0.622 -0.714 -0.252
Texture     0.297 -0.478  0.816
Smoothness  0.721  0.441 -0.115
Compactness 0.949  0.097  0.027
Concavity   0.934 -0.172 -0.039
ConcavePts  0.923 -0.308 -0.159
Symmetry    0.671  0.380  0.077
FracDim     0.466  0.783  0.222

bc_loadings_df <- (data.frame(bc_pca$var$coord) %>%
                     rownames_to_column('Variable'))

bc_base <- (ggplot(bc_loadings_df) +
              theme_classic(base_size = 14) +
              labs(y = '', x = 'Loading') +
              geom_vline(aes(xintercept = 0), linetype = 2))

l1 <- (bc_base +
         geom_col(aes(y = reorder(Variable, Dim.1), x = Dim.1)) +
         labs(title = 'Dimension 1'))
l2 <- (bc_base +
         geom_col(aes(y = reorder(Variable, Dim.2), x = Dim.2)) +
         labs(title = 'Dimension 2'))
l3 <- (bc_base +
         geom_col(aes(y = reorder(Variable, Dim.3), x = Dim.3)) +
         labs(title = 'Dimension 3'))

c1 <- fviz_contrib(bc_pca, choice = 'var', axes = 1)
c2 <- fviz_contrib(bc_pca, choice = 'var', axes = 2)
c3 <- fviz_contrib(bc_pca, choice = 'var', axes = 3)

(c1 + c2 + c3) / (l1 + l2 + l3)

PC1 (53.6%): PC1 mostly measures how irregular the cell shape is. Compactness, Concavity, and ConcavePts are the biggest contributors. Higher scores mean more irregular, indented shapes, while lower scores mean smoother, rounder cells. This could be the main component separating malignant and benign cells.

PC2 (22.8%): PC2 mainly compares FracDim and Radius. It separates smaller cells with more complex borders from larger cells with simpler borders. This captures information that PC1 doesn’t explain very well.

PC3 (10.3%): PC3 is mostly driven by Texture, which contributes about 81%. It measures variation in gray-scale values, with higher scores meaning more texture variation. This makes PC3 worth keeping since it captures information the first two PCs miss.

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

bc_update <- bind_cols(bc_cells, bc_pca$ind$coord[, 1:3])

(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_size_continuous(range = c(0.5, 4)) +
    labs(color = 'Diagnosis', size = 'PC3 score',
         title = 'Biplot of first two PCs'))

(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) +
    theme_classic(base_size = 14))

• Malignant cells have much higher PC1 scores than benign cells. This suggests malignant cells have more irregular, indented shapes, while benign cells tend to be smoother and rounder. Malignant cells also have lower PC2 scores, suggesting they tend to be larger.

• PC1 separates the groups best by far. The boxes don’t overlap, showing a clear difference between malignant and benign cells. PC2 provides some separation, but there is considerable overlap. PC3 does very little to distinguish the groups.