/
Vig_04_Scores_Loadings.Rmd
190 lines (152 loc) · 14.7 KB
/
Vig_04_Scores_Loadings.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
---
title: "Understanding Scores and Loadings"
author:
- David T. Harvey^[Professor of Chemistry & Biochemistry, DePauw University, Greencastle IN USA., harvey@depauw.edu]
- Bryan A. Hanson^[Professor Emeritus of Chemistry & Biochemistry, DePauw University, Greencastle IN USA., hanson@depauw.edu]
date: "`r Sys.Date()`"
output:
bookdown::html_document2: # use for pkgdown site
# bookdown::pdf_document2: # use for CRAN to keep size down; breaks GHA
toc: true
toc_depth: 2
fig_caption: true
number_sections: false
mathjax: default
css: vignette.css
vignette: >
%\VignetteIndexEntry{Vignette 04: Understanding Scores and Loadings}
%\VignetteKeywords{PCA}
%\VignettePackage{LearnPCA}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
#link-citations: true
#bibliography: PCA.bib
#biblio-style: plain
pkgdown:
as_is: true
---
```{r SetUp, echo = FALSE, eval = TRUE, results = "hide"}
# R options & configuration:
set.seed(13)
suppressPackageStartupMessages(library("knitr"))
suppressPackageStartupMessages(library("plotrix"))
suppressPackageStartupMessages(library("latex2exp"))
# source functions
source("scores_loadings_functions.R")
# colors and color names
axis2_col = "#3db7ed"
axis2_colname = "light blue"
axis1_col = "#f748a5"
axis1_colname = "pink"
loadings_col = "#359b73"
loadings_colname = "green"
origdata_col = "#000000"
origdata_colname = "black"
# Stuff specifically for knitr:
opts_chunk$set(eval = TRUE, echo = FALSE, results = "hide")
opts_knit$set(eval.after = "fig.cap")
```
<!-- This chunk inserts common info about all the vignettes -->
```{r, echo = FALSE, results = "asis"}
res <- knitr::knit_child("top_matter.qmd", quiet = TRUE)
cat(res, sep = '\n')
```
```{r prepData}
# create data sets for use in vignette
data_set_1 = generate_data(size = 10)
rotate0 = rot_axes(angle = 0, file = data_set_1)
rotate15 = rot_axes(angle = 15, file = data_set_1)
rotate20 = rot_axes(angle = 20, file = data_set_1)
rotate30 = rot_axes(angle = 30, file = data_set_1)
rotate40 = rot_axes(angle = 40, file = data_set_1)
rotate45 = rot_axes(angle = 45, file = data_set_1)
rotate60 = rot_axes(angle = 60, file = data_set_1)
rotate75 = rot_axes(angle = 75, file = data_set_1)
rotate80 = rot_axes(angle = 80, file = data_set_1)
```
# Introduction
In the vignette A [Conceptual Introduction to PCA](#top-matter), we used a small data set---the relative concentrations of 13 elements in 180 archaeological glass artifacts---to highlight some key features of a principal component analysis. We learned, for example, that a PCA analysis could reduce the complexity of the data from 13 variables to three, which we called the principal components. We also explored how we can use the scores returned by a PCA analysis to assign each of the 180 samples into one of four groups based on the first two principal components, and we learned how we can use the loadings from a PCA analysis to identify the relative importance of the 13 variables to each of the principal components. The vignettes [The Math Behind PCA](#top-matter) and [PCA Functions](#top-matter) explained how we extract scores and loadings from the original data and introduced the various functions within `R` that we can use to carry out a PCA analysis. None of these vignettes, however, explain the relationship between the original data and the scores and loadings we extract from that data by a PCA analysis. As we did in the vignette [Visualizing PCA in 3D](#top-matter), we will use visualizations to help us understand the origin of scores and loadings.
# A Small Data Set
For this vignette we will use a small data set that consists of eight samples and two variables. We are limiting ourselves to two variables so that we can plot the data in two dimensions and we are limiting ourselves to eight samples so that our plots are reasonably uncluttered and easy to view. Figure \@ref(fig:fig-origData) shows a plot of the original data and Table \@ref(tab:table-origData) shows the individual values for each sample and variable.
```{r fig-origData, fig.align="center", fig.cap="Scatterplot of variable 1 and variable 2 for our eight samples. The dashed lines are the original axes with variable 1 shown in pink and variable 2 shown in blue."}
plot_rot_axes(rotate0, show_rotated = FALSE, show_full_legend = FALSE,
show_scores = FALSE, show_simple_legend = FALSE,
show_loadings = FALSE, show_title = FALSE, show_projections = FALSE)
legend(x = "topleft", legend = c("variable 1", "variable 2"),
lty = 2, lwd = 1, col = c(axis1_col,axis2_col), bty = "n")
```
```{r table-origData, results = "asis"}
variable2 = data_set_1$xrot
variable1 = data_set_1$yrot
samples = seq(1,8)
df = data.frame(samples,variable1,variable2)
colnames(df) = c("sample", "variable 1", "variable 2")
kable(df, caption = "Individual values for the variables and samples.", digits = 2)
var2 = var(variable2)
var1 = var(variable1)
total_variance = var1 + var2
percent_var2 = 100*var2/total_variance
percent_var1 = 100*var1/total_variance
```
A cursory examination of Figure \@ref(fig:fig-origData) shows us that there is a positive correlation between the two variables; that is, a general increase in the value for one variable results in a general increase in the value for the other variable. As you might expect from other vignettes, this correlation between the two variables suggests that a single principal component is likely sufficient to explain the data.
Another way to examine our data is to consider the relative dispersion in the values for each variable. An examination of the data in Table \@ref(tab:table-origData) shows that variable 2 spans a greater range of values (from `r round(min(variable2), digits = 2)` to `r round(max(variable2), digits = 2)`) than do the values for variable 1 (from `r round(min(variable1), digits = 2)` to `r round(max(variable1), digits = 2)`). We can treat this quantitatively by reporting the variance[^1] for each variable and the total variance, which is the sum of the variances for the two variables:
- variable 1's variance: `r round(var1, digits = 2)`
- variable 2's variance: `r round(var2, digits = 2)`
- total variance: `r round(total_variance, digits = 2)`
Variable 2 accounts for `r round(percent_var2, digits = 1)`% of the total variance and variable 1 accounts for the remaining `r round(percent_var1, digits = 1)`% of the total variance. It is important to note that these values refer to variance measured relative to the original axis system. If we change the axis system (as we will do momentarily), *the variance relative to each new axis will be different because the position of the points relative to the axes are different.*
# Rotating the Axes
As outlined in the vignette [Visualizing PCA in 3D](#top-matter), a principal component analysis essentially is a process of rotating our original set of $n$ axes, which correspond to the $n$ variables we measured, until we find a new axis that explains as much of the total variance as possible. This becomes the first principal component axis. We then project the data onto the $n - 1$ dimensional surface that is perpendicular to this axis and repeat this process of rotation and projection until the original $n$ axes are replaced with a new set of $n$ principal component axes.
For a system with just two variables, this process amounts to a simple rotation of the two axes. Suppose we arbitrarily rotate the axes in Figure \@ref(fig:fig-origData) clockwise by 20$^\circ$ and project the data points onto the new axes. Figure \@ref(fig:rotate20) shows the result where the dashed lines are the original axes associated with variable 1 and variable 2, and the solid lines are the rotated axes, which we identify for now as axis A and axis B to clarify that they are not specifically associated with one of the two variables, nor with the principal components we eventually hope to find (think of axes A and B as "proposed principal components"). It is easy to see---both visually and quantitatively---that the total variance is more evenly distributed between the two rotated axes than was the case for the original two axes. Even more importantly, the variance associated with axis A has increased from 22.8% to 55.4%. One can see this graphically by looking at the span of the pink and blue projection points on each axis.
```{r rotate20, fig.align="center", fig.cap="Rotating the original axes by 20$^\\circ$ gives the axes shown as solid lines. The small points in pink and in blue are the projections of the original data onto the rotated axes. The gray lines serve as a guide to illustrate the projections for one data point."}
plot_rot_axes(rotate20, show_loadings = FALSE, show_simple_legend = FALSE,
show_full_legend = FALSE, show_full_axes_legend = TRUE,
show_scores = TRUE, show_projections = TRUE)
```
If we continue to rotate the axes, we discover that an angle of rotation of 60$^\circ$ yields an axis (axis A) that explains almost 99% of the total variance in our eight samples. At this point we have maximized the variance for this axis, because further rotation will reduce the variance. As a result, we can now call this axis the first principal component axis. The remaining axis, which is perpendicular to the first principal component axis, is the second principal component axis.[^5]
```{r rotate60, fig.align="center", fig.cap="Rotating the original axes by 60$^\\circ$ maximizes the variance along one of the two rotated axes; the line in pink is the first principal component axis and the line in blue is the second principal component axis. The gray lines serve as a guide to illustrate the projections for one data point."}
plot_rot_axes(rotate60, show_loadings = FALSE, show_simple_legend = FALSE,
show_full_legend = TRUE, show_full_axes_legend = FALSE,
show_scores = TRUE, show_projections = TRUE)
```
If you would like to further explore the interplay of the data, the rotation of the proposed axes and the variance, please try the shiny app by running `PCsearch()`.
# Scores
In Figure \@ref(fig:fig-origData), each or our eight samples appears as a point in space defined by its position along the original axes for variable 1 and variable 2. After completing the rotation of the axes, each of the eight samples now appears as a point in space defined by its position along the recently discovered two principal component axes. These coordinates are the scores returned by the PCA analysis. Table \@ref(tab:table-scores) provides the scores for our eight samples in the columns labeled PC 1 and PC 2; also shown are the values for variable 1 and variable 2.
```{r table-scores, results = "asis"}
variable2 = data_set_1$xrot
variable1 = data_set_1$yrot
samples = c(1:8, "% variance")
pc1 = sqrt(rotate60$axis1_x^2 + rotate60$axis1_y^2) * sign(rotate60$axis1_x)
pc2 = sqrt(rotate60$axis2_x^2 + rotate60$axis2_y^2) * sign(rotate60$axis2_x)
# compute variances for a summary row
v1 <- var(variable1)
v2 <- var(variable2)
vpc1 <- var(pc1)
vpc2 <- var(pc2)
# add variances to data
variable1 <- c(variable1, 100*v1/(v1 + v2))
variable2 <- c(variable2, 100*v2/(v1 + v2))
pc1<- c(pc1, 100*vpc1/(vpc1 + vpc2))
pc2<- c(pc2, 100*vpc2/(vpc1 + vpc2))
# put it all together
df = data.frame(samples,variable1,variable2,pc1, pc2)
colnames(df) = c("sample", "variable 1", "variable 2", "PC 1", "PC 2")
kable(df, caption = "Coordinates for each of the eight samples, in the original axis system and in the principal component axis system.", digits = 2)
```
# Loadings
Although the scores define the location of our samples in the coordinates of the newly discovered principal component axes, they do not tell us where these new axes are located. If we want to reconstruct the original data from the results of a PCA analysis (a process described in more detail in the vignette [Step-by-Step PCA](#top-matter)), we must know both where the samples lie relative to the principal component axes **and** the location of the principal component axes relative to the original axes. As shown in Figure \@ref(fig:rotate60-loadings), each principal component axis is defined by the cosine of its angle of rotation relative to each of the original axes. For example, the angle of rotation of PC 1 to the axis for variable 1, identified here as $\Theta$, is --60$^\circ$, which gives its loading as $\cos(-60) = 0.500$.[^2] The angle of rotation of PC 1 to the axis for variable 2, or $\Phi$, returns a loading of $\cos(30) = 0.866$.[^4]
```{r rotate60-loadings, fig.align="center", fig.cap="Illustration showing the loadings for the first principal component axis."}
plot_rot_axes(rotate60, show_loadings = TRUE, show_simple_legend = FALSE,
show_full_legend = TRUE, show_full_axes_legend = FALSE,
show_scores = TRUE, show_projections = TRUE)
```
Since individual loadings are defined by a cosine function, they are limited to the range --1 (an angle of rotation of 180$^\circ$) to +1 (an angle of rotation of 0$^\circ$).[^3] The sign of the loading indicates how a variable contributes to the principal component. A positive loading indicates that a variable contributes to some degree to the principal component, and a negative loading indicates that its absence contributes to some degree to the principal component. The larger a loading's relative magnitude, the more important is its presence or absence to the principal component.
<!-- Insert refer_to_works_consulted document -->
```{r, echo = FALSE, results = "asis"}
res <- knitr::knit_child("refer_to_works_consulted.md", quiet = TRUE)
cat(res, sep = '\n')
```
[^1]: Variance is defined as $$\frac{\sum({x_i - \bar{x}})^2}{n - 1}$$ for a set of $n$ values of $x$.
[^2]: The optimal angle for this data set is actually 62$^\circ$, but we'll use the rounded value for discussion.
[^3]: And it follows if the principal component axis aligns with an original axis, the loading will be $\cos(0) = 1.00$, and if it is perpendicular, the loading will be $\cos(90) = 0.00$.
[^4]: At this point one should ask which of these two options is the correct one? Should we take our loading to be relative to the variable 1 axis, or the variable 2 axis? It doesn't really matter, as long as once we make the choice the rest of the computations are relative to the same choice. This is related to the fact that the signs of the loadings can change between programs and builds, discussed further in the [PCA Functions](#top-matter) vignette.
[^5]: When the axis is optimally aligned with the data, the maximum variance is explained. Ideally, the spread of data long this axis represents to a high degree the phenomenon under study. Put another way, it represents the signal of interest. You can look at PCA as optimizing the signal to noise ratio along the first principal component axis, with less signal and more noise along the second principal component axis, and so on with each succeeding axis. This is reflected in the scree plot.