Dimension Reduction

MSSC 6250 Statistical Machine Learning

Dr. Cheng-Han Yu
Department of Mathematical and Statistical Sciences
Marquette University

Unsupervised Learning

Unsupervised Learning

• Supervised Learning: response \(Y\) and features \(X_1, X_2, \dots, X_p\) measured on \(n\) observations.

• Unsupervised Learning: only features \(X_1, X_2, \dots, X_p\) measured on \(n\) observations.

  • Not interested in prediction (no response to be predicted)
  • Discover any interesting pattern or relationships among these features.

• Estimate the density, covariance, graph (network), etc. of \(\mathbf{X}\)

• Dimension reduction for effective data visualization or extracting most important information those features contain.
  • plot a bunch of points of \(\boldsymbol{x} = (x_1, x_2, \dots, x_p)\) in a 2-D scatter plot (manifold). (reduce dimension from \(p\) to 2)
  • use 2 variables to explain most variations or represents high data density in the \(p\) variables

• Clustering discovers unknown subgroups/clusters in data
  • find 3 sub-groups of people based on variables income, occupation, age, etc

Background: Dimensions

One-Dimension (1D) Number line

# A tibble: 50 × 1
   English
     <int>
 1      41
 2      65
 3      55
 4      94
 5      66
 6      85
 7      44
 8      44
 9      67
10      73
# ℹ 40 more rows

• Let’s go back to your first grade. You remember that one-dimension equals a number line.
• Now suppose we have 50 English grades. The first 3 scores are 41, 65, and 55.
• We can plot these values on the number line just like we did in the elementary school.
• 1st student with score 41 it’s a dot at 41, 2nd student with score 65, it’s a dot at 65, and so on.

One-Dimension (1D) Number line: Uniform students

# A tibble: 50 × 1
   English
     <int>
 1      41
 2      65
 3      55
 4      94
 5      66
 6      85
 7      44
 8      44
 9      67
10      73
# ℹ 40 more rows

• If you plot all students scores, we might see something like this: an uniform distribution of English scores.

1D Number line: Non-uniform students

# A tibble: 50 × 1
   English
     <int>
 1      77
 2      78
 3      81
 4      78
 5      52
 6      62
 7      47
 8      58
 9      43
10      59
# ℹ 40 more rows

• Or we might get a non-uniform distribution of English scores. Some students obtain a very high grade, but some students are falling behind and get a very low grade.

Two-Dimensions (2D) X-Y Scatter plot: High Correlated

English and Math measure an overall academic performance.

# A tibble: 50 × 2
   English  Math
     <int> <dbl>
 1      41  33.2
 2      65  63.6
 3      55  44.6
 4      94  95  
 5      66  65.6
 6      85  73.1
 7      44  46.6
 8      44  51  
 9      67  69.4
10      73  66.5
# ℹ 40 more rows

• Now let’s go to sixth grade when we learned about 2 dimensional graphs.
• Now we have two axes instead of one, and now we can plot data from two different subjects instead of just one.
• Additional to English scores, here we also have Math scores of some students.
• And we all know how to map a English and math score pair to a point on a 2D plot.
• If we put all students English and Math scores on the plot, we might see something like this that English and Math scores are positively correlated, meaning that a student having high English score tends to have high Math score too.
• This might be due to the fact that English and Math measure an overall academic performance. A good student tends to have both high English and Math scores.

Two-Dimensions (2D) X-Y Scatter plot: No correlated

English and Math measure different abilities.

# A tibble: 50 × 2
   English  Math
     <int> <int>
 1      41    27
 2      65    47
 3      55    32
 4      94    44
 5      66    20
 6      85    30
 7      44    16
 8      44    72
 9      67    86
10      73    82
# ℹ 40 more rows

• Or we might see English and Math scores are not correlated, meaning that a high English score does not tell us anything about Math score, whether it is high or low.
• Because English and Math might measure different abilities.
• English measures verbal and communication skills and Math measures logic and quantitative skills

Three-Dimensions (3D) X-Y-Z Scatter plot

# A tibble: 50 × 3
   English  Math Biology
     <int> <dbl>   <dbl>
 1      41  33.2    39.2
 2      65  63.6    61.6
 3      55  44.6    41.6
 4      94  95      92  
 5      66  65.6    73.6
 6      85  73.1    71.1
 7      44  46.6    56.6
 8      44  51      56  
 9      67  69.4    79.4
10      73  66.5    59.5
# ℹ 40 more rows

• OK. Now maybe in college, we started drawing 3D graph, having height, width and depth.
• With 3 separate axes, we can now plot data from three different subjects.
• So now we have three subjects, English, Math and Biology.
• With the same logic as 2D plot, we put English score on X-axis, Math score on Y-axis, and biology score on Z-axis, then draw lines perpendicular to each axis to put the point in the 3D space.
• And this is our 3D scatter plot.
• So, if we have one subject data, we can make a 1D graph.
• If we have two subject data, we can make a 2D graph.
• If we have three subject data, we can make a 3D graph.

Four-Dimension (4D) X-Y-Z-? Scatter plot

# A tibble: 50 × 4
   English  Math Biology History
     <int> <dbl>   <dbl>   <dbl>
 1      41  33.2    39.2      51
 2      65  63.6    61.6      53
 3      55  44.6    41.6      63
 4      94  95      92        83
 5      66  65.6    73.6      51
 6      85  73.1    71.1      74
 7      44  46.6    56.6      34
 8      44  51      56        33
 9      67  69.4    79.4      76
10      73  66.5    59.5      74
# ℹ 40 more rows

• But what happens if we have data from 4 academic subjects?? here English, Math, Biology and History.
• Well we need 4D graph. The problem is You can’t draw a 4D plot on paper!
• What if a student has 20 different grades from 20 different courses?? We need a 20 dimensional graph?
• But what is that?? There is no way we can draw that.

How about Pair Plots?

• Plotting all 2-D scatter plots of all possible pairs might be a solution to check the relationship between variables, or exploratory data analysis.
• Like here we see basically any two variables or two subject scores are positively correlated.

Tooooo Many Pair Plots!

• If we have \(p\) variables, there are \({p \choose 2} = p(p-1)/2\) pairs.
• If \(p = 10\), we have 45 such scatter plots to look at!
• In practice, we may encounter over 100 variables!!

• But still, there is a problem.
• If you have \(p\) variables, you are gonna have \({p \choose 2} = p(p-1)/2\) pairs.
• If \(p = 10\), you have 45 such scatter plots to look at!
• Not to mention that in real data science work, you may encounter over 100 variables!!

Dimension Reduction

• One variable represents one dimension.

• With many variables in the data, we live in a high dimensional world.

GOAL:

• Find a low-dimensional (usually 2D) representation of the data that captures as much of the information all of those variables provide as possible.

• Use two created variables to represent all \(p\) variables, and make a scatter plot of the two created variables to learn what our observations look like as if they lived in the high dimensional space.

Why and when can we omit dimensions?

• So in order to meaningfully represent the relationship of all variables, we need a technique, Dimension Reduction.
• In mathematics, One variable represents one dimension, so with many variables in the data, we live in a high dimensional world.
• We would like to find a low-dimensional (usually 2D) representation of the data that captures as much of the information all of those variables provide as possible.
• We use two created variables to represent all \(p\) variables, and make a scatter plot of the two created variables to learn what our observations look like as if they lived in the high dimensional space.
• Of course, it’s not always a good idea to just use two variables to represent all \(p\) variables. Some information will be missing when we just use two variables to represent or explain all the relationships of \(p\) variables.
• But sometimes, a low-dimensional representation looks very like a high dimensional space, and does not lose much information. In this situation, a low-dimensional representation is very useful.
• So Let’s see why and when can we omit dimensions?

Variation mostly from One Variable

• Almost all of the variation in the data is from left to right.

• Let’s suppose we have two variables, and their data look like this.
• Here we see almost all of the variation in the data from left to right.
• It means that for variable X, some observations have low values, and some observations have high values.
• But it looks like all observations have variable Y at about the same level.
• Since X varies more than Y, having larger variation, X contains more information than Y in the data.

Variation mostly from One Variable

• If we flattened the data, the graph would not look much different.

• If we flattened the data, removing the up and down variation, our graph would not look much different from what it look like before.
• Most of the variation, or information is still kept in the data.

Variation mostly from One Variable

• If we flattened the data, we could graph it with a 1D number line!

• Now, with up-and-down variation removed, the variable Y becomes redundant, and we can just represent the flattened data using 1D single number line.

Variation mostly from One Variable

• Both graphs say “the important variation is left to right.”

• So in this case, we can display 2D data on a 1D plot without losing too much information!
• Both graphs say “the important variation is left to right.”
• Another example is watching TV. TV is a 2D thing, but we watch TV for 3D shows. The 2D TV represents the 3D shows, programs or drama very well because the 2D representation does not lose much information the 3D shows provide.
• So through the example, we know some dimensions are more important than others. Like here Variable X is more important than variable Y because variable X explains most of the variation stored in the data.

Principal Component Analysis (PCA)

Idea of PCA

• PCA is a dimension reduction tool that finds a low-dimensional representation of a data set that contains as much as possible of variation.

• Each observation lives in a high-dimensional space (lots of variables), but not all of these dimensions (variables) are equally interesting/important.

• The concept of interesting/important is measured by the amount that the observations vary along each dimension.

• PCA is a dimension reduction tool that finds a low-dimensional representation of a data set that contains as much as possible of variation stored in the data set.
• As we’ve seen before, each of the observations lives in a high-dimensional space, meaning that each observation has lots of variables associated with it, but not all of these dimensions (variables) are equally interesting/important.
• The concept of interesting/important is measured by the amount that the observations vary along each dimension.
• A characteristic or attribute of observations is called variable because its value varies from sample to sample.
• If the variable does not vary, it becomes an irrelevant or un-important variable because we cannot use the variable to differentiate or distinguish observations. If everyone in this class gets grade A, then the data science grade is not an important variable to learn which students perform academically better than others.

PCA Illustration: 2 Variable Example

# A tibble: 50 × 2
   English  Math
     <int> <dbl>
 1      41  33.2
 2      65  63.6
 3      55  44.6
 4      94  95  
 5      66  65.6
 6      85  73.1
 7      44  46.6
 8      44  51  
 9      67  69.4
10      73  66.5
11      47  58.9
12      66  66.5
13      57  51  
14      57  33.2
15      77  83.6
16      83  78.8
# ℹ 34 more rows

• We use the English and Math score example to illustrate how PCA works.
• From the data, we have a 2D scatter plot like this.
• With the data, we can have the average of English score and the average of math score, shown as the red squared point in the scatter plot.

Step 1: Shift (or standardize) the Data

• So the two variables have both mean 0. If the variables are measured in a different unit, consider standardization, \(\frac{x_i - \bar{x}}{s_x}\).
• Shifting does not change how the data points are positioned relative to each other.

• PCA first shift or standardize our data so that the center is on top of the origin (0, 0) in the graph.
• If the two variables are measured in a different scale, standardization is highly recommended. Here because both English and Math scores are measured using the same scale. We only do the shifting, and no scaling.
• After standardization, all variables will use the same scale to measure their variation, and that makes more sense because PCA wants to find some axis or coordinate lower dimensional representation that contains the most variation, and we don’t want the variation depends on variable’s measurement units.
• Notice that Shifting does not change how the data points are positioned relative to each other. The scattering shape remains the same.

Step 2: Find a Line that Fits the Data the Best

• Start with a line going through the origin.
• Rotate the line until it fits the data as well as it can, given that it goes through the origin.

• The second step of PCA is going to fit a line to the centered or normalized data set.
• To do this, we start by drawing a random line that goes through the origin.

Step 2: Find a Line that Fits the Data the Best

• Start with a line going through the origin.
• Rotate the line until it fits the data as well as it can, given that it goes through the origin.

• Then we rotate the line until if fits the data as well as it can, given that it has to go through the origin.

Step 2: Find a Line that Fits the Data the Best

• Start with a line going through the origin.
• Rotate the line until it fits the data as well as it can, given that it goes through the origin.

• Ultimately, this line fits best.
• So let’s see what makes this line the best fit, and the meaning of this line.

The Meaning of the Best line

• Principal Component 1 (PC1): maximizes the variance of the projected points.

• PC1 is the line in the Eng-Math space that is closest to the \(n\) observations

  • PC1 minimizes the sum of squared distances between the data points and the PC1.

• PC1 is the best 1D representation of the 2D data

• To quantify how good this line fits the data, PCA projects the data onto it.
• Then the idea is that we can either measure the distances from the data to the line and try to find the line that minimizes those distances.
• Or we can try to find the line that maximizes the distances from the projected points to the origin.
• The two criteria are equivalent.
• DEMO
• So we are find the line that maximizes the variation of the projected points from the data points onto the line!
• The best line is called Principal Component 1 (PC1), which maximizes the sum of squared distances between projected points and the origin.
• Regression line: minimizes the sum of squared residuals (vertical lines from the data points to the line)
• Principal Component 1 (PC1): maximizes the sum of squared distances between between projected points and the origin.
• PC1 is the line in the Eng-Math 2 dimensional space that is closest to the \(n\) observations, i.e., PC1 minimizes the sum of squared distances between the data points and the PC1.
• PC1 is the best 1D representation of the 2D data

The Meaning of the Best line

Source: https://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues

• I am going to conclude the idea of PCA using this gif.
• This gif shows you the idea of finding principal components.
• As the line rotates, you can see the locations of the projected points on the line keeps changing as well.
• And the PC1 is the line that maximizes the variation of the projected points.
• Also, the PC1 is the line that minimizes the distance between the data points and the line. In the figure, the sum of those red lines will be the smallest.
• Questions?

PC1 and PC2

• The data points are also spread out a little above and below the PC1.
• There are some variation that is not explained by the PC1.
• Find the second PC, PC2, that
  • explains the remaining variation
  • is the line through the origin and perpendicular to PC1.

• Along with the PC1, the data points are also spread out a little above and below the PC1.
• There are some variation of the two variables that is not explained by the PC1.
• To find another PC that explains the remaining variation, we find the second PC, called PC2 that is the line through the origin and perpendicular to PC1.

Linear Combinations

• PC1 = 0.68 \(\times\) English \(+\) 0.74 \(\times\) Math
• PC2 = 0.74 \(\times\) English \(-\) 0.68 \(\times\) Math
• PC1 is like an overall intelligence index as it is a weighted average combining verbal and quantitative abilities.
• PC2 accounts for individual difference in English and Math scores.
• \(0.68^2 + 0.74^2 = 1\) (Pythagorean theorem)
• The combination weights 0.68, 0.74, etc are called PC loadings.

• Now let’s look at PC1 and PC2 a little more carefully.
• First PC1 and PC2 are just a vector in 2 dimensional space, right?
• In other words, they are linear combinations of two standard basis, here our English axis and Math axis.
• PCA help us get the linear combinations.
• Here PC1 = 0.68 \(\times\) English + 0.74 \(\times\) Math
• PC2 = 0.74 \(\times\) English - 0.68 \(\times\) Math
• So to make PC1, we mix 0.68 part of English score with 0.74 parts of Math score.
• One unit of PC1 consists of 0.68 parts of English and 0.74 parts of Math.
• And because the weight of math score is a little bit larger, Math score is a little bit more important when it comes to describing how the data are spread out.
• PC1 can be viewed as an overall intelligence index because it is an weighted average combining both verbal and quantitative reasoning abilities.
• PC2 accounts for individual difference in English and Math scores.
• \(0.68^2 + 0.74^2 = 1\) (Pythagorean theorem)
• The combination weights 0.68, 0.74, etc are called loadings of PC.

Rotate Everything so that PC1 is Horizontal

1D representation

• PC1 is our 1D number line that explains the most variation contained in 2D data using a 1D line.
• Points on the PC1 are the projected points of data onto PC1.

2D representation

• The new coordinates PC1 and PC2 are ordered by variation size of the English and Math scores

• If we make PC1 and PC2 the new X-Y coordinates, the new coordinates PC1 and PC2 are ordered by variation size of the English and Math scores.
• PC1 is our 1D number line that explains the most variation contained in 2D data using a 1D line.
• Points on the PC1 are the projected points of data onto PC1.
• PC1 explains the most variation, and PC2 explains the second most variation.
• And PC2 explains the remaining variations that are not explained by PC1.
• Here, because we only have 2 dimensions, we have at most 2 PCs.

Variation

• If the variation for PC1 is \(17\) and the variation for PC2 is \(2\), the total variation presented in the data is \(17+2=19\).

• PC1 accounts for \(17/19 = 89\%\) of the total variation, and PC2 accounts for \(2/19 = 11\%\) of the total variation.

• OK how do we quantify variation. Let’s give it a definition.
• Variation for PC1 \(= \frac{\text{Sum of squared distances of projected points on PC1}}{n-1}\)
• Variation for PC2 \(= \frac{\text{Sum of squared distances of projected points on PC2}}{n-1}\)
• If the variation for PC1 is \(17\) and the variation for PC2 is \(2\), the total variation presented in the data is \(17+2=19\).
• PC1 accounts for \(17/19 = 89\%\) of the total variation, and PC2 accounts for \(2/19 = 11\%\) of the total variation.

How about 3 or More Variables?

• PC1 spans the direction of the most variation

• PC2 spans the direction of the 2nd most variation

• PC3 spans the direction of the 3rd most variation

• PC4 spans the direction of the 4th most variation

• If we have \(n\) observations and \(p\) variables (dimensions), there are at most \(\min(n - 1, p)\) PCs.

• OK. We use a 2D data to illustrate PCA. How about 3 or More Variables?
• Well the idea of 2D data can be exactly applied to 3 or more dimensional data.

US Arrest Data in 1973

dim(USArrests)

[1] 50  4

head(USArrests, 16)

            Murder Assault UrbanPop Rape
Alabama       13.2     236       58   21
Alaska        10.0     263       48   44
Arizona        8.1     294       80   31
Arkansas       8.8     190       50   20
California     9.0     276       91   41
Colorado       7.9     204       78   39
Connecticut    3.3     110       77   11
Delaware       5.9     238       72   16
Florida       15.4     335       80   32
Georgia       17.4     211       60   26
Hawaii         5.3      46       83   20
Idaho          2.6     120       54   14
Illinois      10.4     249       83   24
Indiana        7.2     113       65   21
Iowa           2.2      56       57   11
Kansas         6.0     115       66   18

• OK. Time to learn how to perform PCA in R.
• The data set used for PCA is USArrests data set in 1973.
• Each observation or subject is a state, and we have 4 features or dimensions, Murder, Assault, UrbanPop, and Rape.

PC Loading Vectors on USArrests

pca_output <- prcomp(USArrests, scale = TRUE)

## rotation matrix provides PC loadings
(pca_output$rotation <- -pca_output$rotation)

          PC1   PC2   PC3    PC4
Murder   0.54  0.42 -0.34 -0.649
Assault  0.58  0.19 -0.27  0.743
UrbanPop 0.28 -0.87 -0.38 -0.134
Rape     0.54 -0.17  0.82 -0.089

• PCs are unique up to a sign change, so -pca_output$rotation gives us the same PCs as pca_output$rotation does. The sign just change the direction, not the angle.

\(\text{PC1} = 0.54 \times \text{Murder} + 0.58 \times \text{Assault} + 0.28 \times \text{UrbanPop} + 0.54 \times \text{Rape}\)

\(\text{PC2} = 0.42 \times \text{Murder} + 0.19 \times \text{Assault} - 0.87 \times \text{UrbanPop} - 0.17 \times \text{Rape}\)

• We have 4 PCs because \(\min(n-1, p) = \min(50-1, 4) = 4\).

• To perform PCA in R, it cannot be easier.
• We just need to use the function prcomp(), and put the data set in the function. Then we get everything we want.
• Here I choose to scale the data because variables are not measured in the same scale. For example, Murder rate and UrbanPop are measured in different units.
• This makes sure that every variable has variance 1, and our analysis is not affected by units.
• The PCA results are stored as a list in pca_output object.
• OK first we can look at the rotation matrix because it provides PC loadings, and so we know what PC1 and PC2 are.
• Changing signs for easier interpretation of PCs
• Those PC loadings define how we rotates the coordinates to obtain the PCs. ???
• Again, PC1 is just a linear combination or weighted average of the 4 variables, same as PC2. ???
• Again, PC1 is just a linear combination or weighted average of the 4 variables, same as PC2.
• PCs are unique up to a sign change, so -pca_output$rotation gives us the same PCs as pca_output$rotation does.
• We have 4 PCs because \(\min(n-1, k) = \min(50-1, 4) = 4\).

PC Scores

• The value of the rotated data, the data values of each PC are stored in pca_output$x

head(pca_output$x <- -pca_output$x, 16) |> round(2)

              PC1   PC2   PC3   PC4
Alabama      0.98  1.12 -0.44 -0.15
Alaska       1.93  1.06  2.02  0.43
Arizona      1.75 -0.74  0.05  0.83
Arkansas    -0.14  1.11  0.11  0.18
California   2.50 -1.53  0.59  0.34
Colorado     1.50 -0.98  1.08  0.00
Connecticut -1.34 -1.08 -0.64  0.12
Delaware     0.05 -0.32 -0.71  0.87
Florida      2.98  0.04 -0.57  0.10
Georgia      1.62  1.27 -0.34 -1.07
Hawaii      -0.90 -1.55  0.05 -0.89
Idaho       -1.62  0.21  0.26  0.49
Illinois     1.37 -0.67 -0.67  0.12
Indiana     -0.50 -0.15  0.23 -0.42
Iowa        -2.23 -0.10  0.16 -0.02
Kansas      -0.79 -0.27  0.03 -0.20

• The value of the rotated data, the data values of each PC are stored in pca_output$x.
• In other words, PC1 column here shows the projected values of observations onto PC1. PC2 column shows the projected values of observations onto PC2, and so on.

Interpretation of PCs

pca_output$rotation

          PC1   PC2   PC3    PC4
Murder   0.54  0.42 -0.34 -0.649
Assault  0.58  0.19 -0.27  0.743
UrbanPop 0.28 -0.87 -0.38 -0.134
Rape     0.54 -0.17  0.82 -0.089

• PCs are less interpretable than original features.
• The first loading vector places approximately equal weight on Assualt, Murder and Rape, with much less weights on UrbanPop.
• PC1 roughly corresponds to a overall serious crime rate.

• The second loading vector places most of its weight on UrbanPop, and much less weight on the other 3 features.
• PC2 roughly corresponds to the level of urbanization.

• Intepretability decreases with the order of PCs.
• So it’s easier to give PC1 a meaningful name than PC2, and PC2 is more meaningful than PC3, and so on. Because the PCs after the first 2 PCs usually explain quite small variation in the data, and some of them may be just noises.
• Let’s see if we can interpret these PCs.
• First keep in mind that PCs are less interpretable than original features. Sometimes we even don’t know how to interpret it, especially for PCs that explain small variations. So this is the price we pay for dimension reduction.
• But let’s look at this example.
• The first loading vector places approximately equal weight on Assualt, Murder and Rape, with much less weights on UrbanPop.
• So PC1 roughly corresponds to a overall serious crime rate because PC1 explains the variations of data caused by those crime variables Assualt, Murder and Rape.
• On the contrary, the second loading vector places most of its weight on UrbanPop, and much less weight on the other 3 features.
• So we can say PC2 roughly corresponds to the level of urbanization.
• So you get the idea, Assualt, Murder and Rape are similar each other because they all are measures of crime rate.
• So when reducing dimensions, we sort of combine the three similar variables together to become a one single index that measures an overall crime rate.
• Urban population measures a totally different thing. So the variation created by this variable cannot be explained well by the crime rate, and it should be absorbed in PC2.

2D Representation of the 4D data

            PC1   PC2   PC3   PC4
Wisconsin -2.06 -0.61 -0.14 -0.18
Wyoming   -0.62  0.32 -0.24  0.16

• Higher value of PC1 means higher crime rates (roughly).

• Higher value of PC2 means higher level of urbanization (roughly).

• And we can show our 2D Representation of the 4D data.
• Higher value of PC1 means higher crime rates (roughly).
• Higher value of PC2 means higher level of urbanization (roughly).
• We may be able to further analysis on the reduced dimension data, for example, we may want to partition the data into 2 clusters.
• One cluster has high crime rate and low urbanization, the other group shows low crime rate and high urbanization.
• One cluster may be the low economically-developed states, and the other high economically-developed states.
• So you see, we can tell lots of stories by analyzing our data.

2D Representation of the 4D data: biplot

biplot(pca_output, xlabs = state.abb, 
       scale = 0)

• Top axis: PC1 loadings
• Right axis: PC2 loadings
• Red arrows: PC1 and PC2 loading vector, e.g., (0.28, -0.87) for UrbanPop.
• Crime-related variables (Assualt, Murder and Rape) are located close to each other.
• UrbanPop is far from the other three.
• Assualt, Murder and Rape are more correlated, and UrbanPop is less correlated with the other three.

• We can simply use the function biplot() to show the 2D Representation of the data.
• This function also provide loading vector of PC1 and PC2 that gives us an idea of which direction means a large value of a variable/feature.
• This is why it is called biplot because we plot two things in one single plot.
• Here, Top axis is for PC1 loadings
• Right axis is PC2 loadings
• Red arrows: PC1 and PC2 loading vector, e.g., (0.28, 0.87) for UrbanPop.
• So NJ and Ca have pretty high urban pop rate because they are large in the UrbanPop arrow direction.
• Crime-related variables (Assualt, Murder and Rape) are located close to each other.
• UrbanPop is far from the other three.
• Assualt, Murder and Rape are more correlated, and UrbanPop is less correlated with the other three.
• Assualt, Murder and Rape sort of point to the same direction as PC1 and UrbanPop points to the same direction as PC2.

Proportion of Variance Explained

(pc_var <- pca_output$sdev ^ 2)

[1] 2.48 0.99 0.36 0.17

(pc_var_prop <- pc_var / sum(pc_var))

[1] 0.620 0.247 0.089 0.043

• PC1 explains \(62\%\) of the variations in the data, and PC2 explains \(24.7\%\) of the variance.
• PC1 and PC2 explain about \(87\%\) of the variance, and the last two PCs explain only \(13\%\).
• 2D plot provides pretty accurate summary of the data.

• Finally I want to talk a little bit about Proportion of Variance Explained.
• In the pca_output, we have SD of each PC.
• We square it to get variance.
• Then if we divided by the sum of variance, we get the Proportion of Variance Explained by each PC.
• So PC1 explains \(62\%\) of the variations in the data, and PC2 explains \(24.7\%\) of the variance.
• PC1 and PC2 explain about \(87\%\) of the variance, and the last two PCs explain only \(13\%\).
• 2D plot provides pretty accurate summary of the data.

Scree Plot

Look for a point at which the proportion of variance explained by each subsequent PC drops off.

Mathematics of PCA

Principal Components

• The \(k\)th PC is

\[Z_k = \phi_{1k}X_1 + \phi_{2k}X_2 + \dots + \phi_{pk}X_p,\] where \(\sum_{j=1}^p\phi_{jk}^2=1\).

• \((\phi_{1k}, \phi_{2k}, \dots, \phi_{pk})'\) is the PC loading vector.

• The PC1 loading vector solves

\[\max_{\phi_{11}, \phi_{21}, \dots, \phi_{p1}} \left\{ \frac{1}{n}\sum_{i=1}^n z_{i1}^2\right\} = \left\{ \frac{1}{n}\sum_{i=1}^n \left( \sum_{j=1}^p \phi_{j1} x_{ij}\right)^2\right\} \quad \text{s.t.} \quad \sum_{j=1}^p\phi_{j1}^2 = 1\]

• Maximize the sample variance of the projected points, or the scores \(z_{11}, z_{21}, \dots, z_{n1}\).

• The PC loading vector defines a direction in feature space along which the data vary the most.

Principal Components

For \(k\)th PC, \(k > 1\),

\[\max_{\phi_{1k}, \phi_{2k}, \dots, \phi_{pk}} \left\{ \frac{1}{n}\sum_{i=1}^n z_{ik}^2\right\} = \left\{ \frac{1}{n}\sum_{i=1}^n \left( \sum_{j=1}^p \phi_{jk} x_{ij}\right)^2\right\} \quad \text{s.t.} \quad \sum_{j=1}^p\phi_{jk}^2 = 1, \text{ and } {\mathbf{z}_m}'\mathbf{z}_k = 0, \, m = 1, \dots, k-1\] where

\(\mathbf{z}_l = (z_{1l}, z_{2l}, \dots, z_{nl})'\)

Low-Rank Approximation

• PCs provide low-dimensional planes that are closest to the observations.

• \(x_{ij} \approx \sum_{m=1}^M z_{im}\phi_{jm}\) with equality when \(M = \min(n-1, p)\)

\[(z_{im}, \phi_{jm}) = \mathop{\mathrm{arg\,min}}_{a_{im}, b_{jm}} \left\{ \sum_{j=1}^p\sum_{i=1}^n\left( x_{ij} - \sum_{m=1}^Ma_{im}b_{jm}\right)^2\right\}\]

Source: ISL Fig 12.2

Scaling the Variables

apply(USArrests, 2, var)

  Murder  Assault UrbanPop     Rape 
      19     6945      210       88 

• If we perform PCA on the unscaled variables, PC1 loading vector will have a large loading for Assault.

• When all the variables are of the same type, no need to scale the variables.

PCA and SVD

PCA is equivalent to singular value decomposition (SVD) of \(\mathbf{X}\)

\[\mathbf{X}_{n\times p} = \mathbf{U}_{n \times n} \mathbf{D}_{n \times p} \mathbf{V}'_{p \times p}\] where \(\mathbf{U}\) and \(\mathbf{V}\) are orthogonal matrices, and \(\mathbf{D}\) is diagonal with diagonal elements singular values \(d_1 \ge d_2 \ge \cdots \ge d_p\).

PC Scores and SVD

• \({\bf Z}_{n \times p} = [{\bf z}_1 \, {\bf z}_2 \, \cdots \, {\bf z}_p] = \mathbf{X}\mathbf{V}\)

• The \(j\)th PC is the \(j\)th column of \({\bf Z}\) given by \({\bf z}_j = (z_{1j}, z_{2j}, \dots, z_{nj})' = {\bf Xv}_j\).

• Project \(\mathbf{X}\) onto the space spanned by \(\mathbf{v}_j\)s

• \(\mathbf{v}_j\)s are loading vectors.

• \({\bf Z}_{n \times p} = \mathbf{U}\mathbf{D}\)

• \({\bf z}_j = d_j\mathbf{u}_j\).

• \(\mathbf{u}_j\)s are the unit PC vectors, and \(d_j\) controls the variation along the \(\mathbf{u}_j\) direction.

Low-Rank Approximation and SVD

\({\mathbf{X}} = {\mathbf{U}}{\mathbf{D}}{\mathbf{V}}'\)

\(\mathbf{A}= \mathbf{U}_{n\times M}\mathbf{D}_{M \times M}\) and \(\mathbf{B}' = \mathbf{V}_{M\times p}'\) are the minimizer of

\[ \min_{\mathbf{A}\in \mathbf{R}^{n \times M}, \mathbf{B}\in \mathbf{R}^{p \times M}} \|\mathbf{X}- \mathbf{A}\mathbf{B}' \|\] where

• \(\mathbf{U}_{n\times M}\) is \(\mathbf{U}\) with the first \(M\) columns
• \(\mathbf{D}_{M \times M}\) is \(\mathbf{D}\) of the first \(M\) rows and columns
• \(\mathbf{V}_{M\times p}'\) is \(\mathbf{V}'\) with the first \(M\) rows

• \(x_{ij} = \sum_{m = 1}^p d_mu_{im}v_{jm}\)

• \(x_{ij} \approx \sum_{m = 1}^M d_mu_{im}v_{jm}\)

PCA and Eigendecomposition

PCA is equivalent to eigendecomposition of \(\mathbf{X}'\mathbf{X}\) or \(\boldsymbol \Sigma= \text{Cov}(\mathbf{X}) = \dfrac{\mathbf{X}'\mathbf{X}}{n-1}\)¹, the covariance matrix of \(\mathbf{X}\).

\[\mathbf{X}'\mathbf{X}= \mathbf{V}\mathbf{D}^2\mathbf{V}' = d_1^2\mathbf{v}_1\mathbf{v}_1' + \dots + d_p^2\mathbf{v}_p\mathbf{v}_p'\]

• Total variation: \(\sum_{j=1}^p \text{Var}(\mathbf{x}_j) = \frac{1}{n-1}\sum_{j=1}^pd_j^2 = p\)

• Variation of \(m\)th PC: \(\text{Var}(\mathbf{z}_m) = \frac{d_m^2}{n-1}\)

1. Unit normal scaling (normalization or standardization) is such that the variance is one, and unit length scaling makes the length one. If \(\mathbf{W}\) is the unit length scaled version, \(\text{Corr}(\mathbf{W}) = \mathbf{W}'\mathbf{W}= \frac{\mathbf{X}'\mathbf{X}}{n-1}\)

Transformed OLS Regression

Transform \({\bf y = X\boldsymbol \beta+ \boldsymbol \epsilon}\) into \[{\bf y = XVV'\boldsymbol \beta+ \boldsymbol \epsilon= Z\boldsymbol \alpha} + \boldsymbol \epsilon\] where \(\mathbf{Z}= \mathbf{X}\mathbf{V}\) and \(\boldsymbol \alpha= \mathbf{V}'\boldsymbol \beta\), or \(\boldsymbol \beta= \mathbf{V}\boldsymbol \alpha\).

• The least- squares estimator \(\hat{\boldsymbol \alpha} = (\mathbf{Z}'\mathbf{Z})^{-1}\mathbf{Z}'\mathbf{y}= \mathbf{D}^{-2}\mathbf{Z}'\mathbf{y}\)

• \(\text{Var}\left(\hat{\boldsymbol \alpha} \right) = \sigma^2 (\mathbf{Z}'\mathbf{Z})^{-1} = \sigma^2 \mathbf{D}^{-2}\)

• A small \(d_j\) means that the variance of \(\alpha_j\) will be large.

• The PC regression combats multicollinearity by using less PCs \((m \ll p)\) in the model.

  • The PCs corresponding to tiny \(d_j\)s (huge variance) are removed and least-squares is applied to the remaining PCs.

• Even though the new variables are orthogonal, the same magnitude of variance due to the ill-conditioning in \(X'X\) is retained. The total variance has merely been redistributed.
• \({\bf y = X\boldsymbol \beta+ \boldsymbol \epsilon}\)
• \({\bf y = XVV'\boldsymbol \beta+ \boldsymbol \epsilon= Z\boldsymbol \alpha} + \boldsymbol \epsilon\) where \({\bf Z =XV}\) and \({\bf \boldsymbol \alpha= V'\boldsymbol \beta}\).
• The least- squares estimator \(\hat{\boldsymbol \alpha} = {\bf (Z'Z)^{-1}Z'y = \Lambda^{-1}Z'y}\)
• \(\var\left(\hat{\boldsymbol \alpha} \right) = \sigma^2 {\bf (Z'Z)}^{-1} = \sigma^2 {\bf \Lambda}^{-1}\)
• A small \(\lambda_j\) means that the variance of \(\alpha_j\) will be large.
• The PC regression combats multicollinearity by using less PCs in the model.
• Assume that the regressors are arranged in order of decreasing eigenvalues, \(\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_k > 0\).
• The PCs corresponding to near-zero eigenvalues are removed and least squares applied to the remaining PCs.
• If all PCs are in the model, it is just a rotation of the regressors. The same magnitude of variance is retained.
• The small variation in the \(z_2\) direction is responsible for the large variance in \(\hat{\alpha}_2\).
• A large \(\lambda_1\) allows variance of \(\hat{\alpha}_1\) to be relatively unaffected by the dependency of \(x_1\) and \(x_2\).

Principal Component Regression (PCR)

• Keep \(m < p\) PCs with the largest \(m\) singular values \(d_1, d_2, \dots, d_m\).
• The least- squares estimator is \[\hat{\boldsymbol \alpha}_m = ({\bf Z} _m ' {\bf Z} _m) ^ {-1} {\bf Z}_m ' {\bf y} = \mathbf{D}_m ^ {-2} {\bf Z}_m '{\bf y}\] where \({\bf Z}_m = [{\bf z}_1 \, {\bf z}_2 \, \cdots \, {\bf z}_m]\) and \({\bf \mathbf{D}}_m = \text{diag}(d_1, d_2, \dots, d_m)\)

• \({\bf b}_{pc} = {\bf V}_m\hat{\boldsymbol \alpha}_m = [{\bf v}_1 \, {\bf v}_2 \, \cdots \, {\bf v}_m] \begin{bmatrix} \hat{\alpha}_1 \\ \hat{\alpha}_2 \\ \vdots \\ \hat{\alpha}_m \end{bmatrix}\)
• \({\bf b}_{pc}\) is biased but has smaller variance than \({\bf b}\).

• Keep the first \(m < k\) PCs.
• The least- squares estimator is \[\hat{\boldsymbol \alpha}_m = ({\bf Z} _m ' {\bf Z} _m) ^ {-1} {\bf Z}_m ' {\bf y} = {\bf \Lambda} _m ^ {-1} {\bf Z}_m '{\bf y}\] where \({\bf Z}_m = [{\bf z}_1 \, {\bf z}_2 \, \cdots \, {\bf z}_m]\) and \({\bf \Lambda}_m = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_m)\)
• To transform back to the original coefficients, since \({\bf \boldsymbol \alpha= V'\boldsymbol \beta}\), \({\bf \boldsymbol \beta= V\boldsymbol \alpha}\)
• \({\bf b}_{pc} = {\bf V}_m\hat{\boldsymbol \alpha}_m = [{\bf v}_1 \, {\bf v}_2 \, \cdots \, {\bf v}_m] \begin{bmatrix} \hat{\alpha}_1 \\ \hat{\alpha}_2 \\ \vdots \\ \hat{\alpha}_m \end{bmatrix}\)
• \({\bf b}_{pc}\) is biased but has smaller variance than \({\bf b}\). Deleting PCs does not imply deletion of any of original regressors. The selected \(m\) PCs do contain information provided by all \(k\) original regressors.

Principal Component Regression

• PCR performs well when the directions in which \(X_1, \dots ,X_p\) show the most variation (the first few PCs) are the directions that are associated with \(Y\).

• Deleting PCs does not imply deletion of any of original regressors. The selected \(m\) PCs contain information provided by all \(k\) original regressors. (Not a feature selection method)

• PCR and ridge regression are closely related. Ridge regression is a continuous version of PCR.

pls::pcr()

library(pls)
set.seed(1)
pcr.fit <- pcr(Salary ~ ., data = Hitters, 
               subset = train,
               scale = TRUE, validation = "CV")
validationplot(pcr.fit, val.type = "MSEP")

https://stats.stackexchange.com/questions/371054/what-is-adjusted-cv-or-bias-corrected-cv

pls::pcr()

pcr.pred <- predict(pcr.fit, x[test, ], ncomp = 5)
mean((pcr.pred - y.test)^2)

[1] 142812

pcr.fit <- pcr(y ~ x, scale = TRUE, ncomp = 5)
summary(pcr.fit)

Data:   X dimension: 263 19 
    Y dimension: 263 1
Fit method: svdpc
Number of components considered: 5
TRAINING: % variance explained
   1 comps  2 comps  3 comps  4 comps  5 comps
X    38.31    60.16    70.84    79.03    84.29
y    40.63    41.58    42.17    43.22    44.90

coefficients are in scaled unit, not original units.

Other Dimension Reduction (Latent Variable) Methods

• Kernel Principal Component Analysis https://ml-explained.com/blog/kernel-pca-explained

• Probabilistic PCA

• Factor Analysis

• Autoencoders

• t-SNE