PCA

Key Packages
require(stats, vegan, tidyverse, ggbiplot)

Contents:

 Introduction
 Review of Eigenanalysis
 Theory
 Darlingtonia Dataset
 The Steps of a PCA
 How Many Principal Components?
 PCA Functions in R
 Using PC Scores
 The Biplot: Visualizing a PCA
 Conclusions
 References

Introduction
We are focusing today on Principal Components Analysis (PCA), which is an eigenanalysis-
based approach. We begin, therefore, by reviewing eigenanalysis (for more details on this topic,
refer to the chapter about Matrix Algebra).

Review of Eigenanalysis
Eigenanalysis is a method of identifying a set of linear equations that summarize a symmetric
square matrix. It yields a set of eigenvalues (λ), each of which has an
associated eigenvector (x). The connection between these terms is expressed in Equation A.16
from Gotelli & Ellison (2004):

Ax = λx

In words, this says that a square matrix A can be multiplied by a vector x and yield the same
values as a scalar value λ multiplied by the same vector x. While this may not sound very
helpful, it is the basis for ordination techniques such as PCA. It means that data can be
rotated, reflected, stretched, or compressed in coordinate-space by multiplying the individual
data points by an eigenvector (x).
Another way to describe eigenanalysis is the process of identifying and computing a new
coordinate system for the data. We conventionally represent data graphically such that each axis
is a response variable. Eigenanalysis allows us to create new axes that are combinations of the
response variables.

Eigenanalysis can only be calculated for certain square symmetric matrices (see below). For
example, here’s a simple square symmetric matrix:
A <- matrix(c(3, 1, 1, 0.5), nrow = 2)
An eigenanalysis of this matrix can be conducted using the eigen() function:
E <- eigen(A); E
The resulting object contains the eigenvalues and eigenvectors; they can be extracted for separate
manipulation:
Evalues <- E$values
Evectors <- E$vectors

One property of the matrices to which eigenanalysis can be applied is that the total variance is
equal to the sum of its diagonal values; this is called the trace of the matrix. For matrix A:
sum(diag(A))

Equivalently, the total variance of the matrix is the sum of its eigenvalues:
sum(Evalues)

The eigenvalues (λ in Equation A.16) represent the variance extracted by each axis. There are as
many eigenvalues as there are rows or columns in the matrix (it doesn’t matter which you
consider, as the matrix is symmetric), and the eigenvalues are sorted in order of decreasing size.
In other words, the first axis always accounts for the most variance, the second axis for the next
most variance, etc. They are usually reported as a percentage of the total variance.

The eigenvectors (x in Equation A.16) are associated with the eigenvalue in the same relative
position – i.e., the first eigenvector (first column) is associated with the first eigenvalue. The
sign of the eigenvector may vary among analyses. For example, an element may be positive in
one analysis but negative in another. These analyses are equivalent; they are just mirror images
of one another.

Each combination of an eigenvalue (λ) and its associated eigenvector (x) is a solution for Gotelli
& Ellison’s Equation A.16 (above). For example, verify this equation by comparing the results
using the first eigenvalue and its associated eigenvector:
A %*% Evectors[,1]
Evalues[1] * Evectors[,1]

An important benefit of eigenanalysis is that the axes of the new coordinate system are
uncorrelated with one another. The weights (eigenvectors) are chosen so that the value for a
sample unit on one axis is unrelated to the value for that same sample unit on another axis.
Theory
PCA was developed by Karl Pearson (1901), who used it to determine racial assignment of
individuals based on multiple biometric measurements (Gotelli & Ellison 2004). Harold
Hotelling (1933) developed the mathematics behind PCA, and Goodall (1954) introduced it to
the ecological literature under the term ‘factor analysis’ (which is now used to indicate a
different analytic technique).

It may be helpful to consider a verbal description of what a PCA does. Imagine your data as a
cloud of data points. PCA begins by identifying the ‘best-fitting straight line’ through the cloud.
This line is oriented so that it explains as much of the variability in the data cloud as possible,
and is called the first principal axis or major axis. Now, imagine rotating the data cloud
around this axis until you can identify a second axis that explains as much of the remaining
variation as possible. Since you are holding the first axis fixed, the second axis will be
perpendicular to (orthogonal to) it. Repeat this process for every subsequent axis: rotate the data
cloud around the axes that have already been identified and determine the axis that explains as
much of the remaining variation as possible while remaining orthogonal to all other axes. Note
that although we can perform this mathematically for any number of dimensions, we cannot
visualize it in more than three dimensions!

PCA identifies new ‘synthetic’ variables (principal components; PCs) that are composites or
blends of the original variables. The key properties of a PCA are:

 Each eigenvalue measures the amount of variance explained by that principal component.
 Each eigenvector is the set of weights assigned to the original variables to produce a
principal component.
 If we keep all of the PCs, we have rotated the data without altering the distances among
observations.
 If we focus on the first few PCs, we have reduced the dimensionality of the dataset.
 Since the eigenvalues are always in declining order of importance, focusing on them in
order means that we are focusing on as much variation as possible to explain in that
number of dimensions.
 The principal components are uncorrelated with (orthogonal to) each other. Therefore,
these components may better meet the assumptions of other analytical techniques
(multiple regression, MANOVA, etc.) than the original variables could.
 The effectiveness of PCA is a function of how strongly correlated the original variables
are:
o If the variables are highly correlated, PCA can very effectively express them in a
smaller number of principal components.
o If the variables are uncorrelated, there by definition is no ‘shared’ information
between them and each principal component is largely associated with one of the
variables.

PCA is a very important technique to understand. More details about PCA can be found in
Borcard et al. (2018; Section 5.3) and Legendre & Legendre (2012; Section 9.1).
PCA is appropriate for many types of data (e.g., LiDAR, morphological data). Summerville et
al. (2006) apply this technique to identify suites of correlated traits among moth species, and
then analyze each suite of traits (principal component) individually. However, PCA is widely
acknowledged to be inappropriate for community-level analyses that involve sample × species
data matrices (see ‘Conclusions’ section for details).

Darlingtonia Dataset
We’ll illustrate the steps of a PCA using a dataset containing detailed morphological
measurements of Darlingtonia californica (cobra lily, or California pitcherplant)
plants. Darlingtonia are carnivorous plants with specialized leaves that function as pitchers
which trap insects.

This dataset is taken from Table 12.1 of Gotelli & Ellison (2004). Measurements were made on
87 plants from four sites in southern Oregon. Ten variables were measured, 7 based on
morphology and 3 based on biomass:

Variable Uni Description

 Name ts

height mm Plant height

[Link] mm Mouth opening

[Link] mm Tube diameter

[Link] mm Keel diameter

[Link] Length of one of the two ‘wings’ that make up the

mm
h fishtail appendage

[Link] Length of the other ‘wing’ that makes up the fishtail

mm
h appendage

wingsprea mm Distance between the two wing tips

hoodmass.g g Dry weight of hood

tubemass.g g Dry weight of tube

wingmass.g g Dry weight of fishtail appendage

The data file (Darlingtonia_GE_Table12.[Link]) is available in the ‘Data files’ folder on the
course website. Save it to the ‘data’ sub-folder within the folder containing your SEFS 502
project.
 We begin by loading the data and creating an object that contains the response variables (i.e.,
omitting the columns that code for site and plant number):
darl <- [Link]("data/Darlingtonia_GE_Table12.[Link]",
header=TRUE)
[Link] <- darl %>%
select(-c(site, plant))

The Steps of a PCA

Now, let’s work through a PCA step-by-step. Note that the first two steps can each be done in
two ways.

Step 1: Normalize Data (optional)

Option A: Normalize data by subtracting the mean and dividing by the standard deviation for
each column. This is done so that all variables are expressed on the same scale. The centroid of
the normalized dataset is therefore at the origin.
[Link] <- scale([Link])

If variables are not on the same scale, they will contribute unequally to the results: a variable
with high variance will contribute much more than a variable with low variance.

Option B: Keep data in raw format.

Step 2: Calculate Variance/Covariance Matrix or

Correlation Matrix
Option A: Calculate the variance/covariance matrix (S) based on the scaled data:
S <- cov([Link])

Option B: Convert raw data to correlation matrix (P).

P <- cor([Link])

Look at S and P – they’re identical! (Aside: this is only true if the variables have been
normalized. If you centered but did not rescale them, each variable would be weighted in
proportion to its variance. In my experience this is uncommon.)
S %>% round(2)
height [Link] [Link] [Link] [Link]
height 1.00 0.61 0.24 -0.03 0.28
[Link] 0.61 1.00 -0.05 -0.35 0.57
[Link] 0.24 -0.05 1.00 0.54 -0.09
[Link] -0.03 -0.35 0.54 1.00 -0.33
[Link] 0.28 0.57 -0.09 -0.33 1.00
[Link] 0.29 0.44 0.01 -0.26 0.84
wingsprea 0.12 0.25 0.10 -0.21 0.64
hoodmass.g 0.49 0.76 -0.11 -0.30 0.61
tubemass.g 0.84 0.72 0.10 -0.25 0.44
wingmass.g 0.10 0.24 0.03 -0.12 0.47
[Link] wingsprea hoodmass.g tubemass.g
wingmass.g
height 0.29 0.12 0.49 0.84
0.10
[Link] 0.44 0.25 0.76 0.72
0.24
[Link] 0.01 0.10 -0.11 0.10
0.03
[Link] -0.26 -0.21 -0.30 -0.25 -
0.12
[Link] 0.84 0.64 0.61 0.44
0.47
[Link] 1.00 0.71 0.49 0.36
0.47
wingsprea 0.71 1.00 0.25 0.16
0.34
hoodmass.g 0.49 0.25 1.00 0.76
0.35
tubemass.g 0.36 0.16 0.76 1.00
0.17
wingmass.g 0.47 0.34 0.35 0.17
1.00

Since S and P are identical, we can work with either one in all subsequent steps.

Step 3: Conduct Eigenanalysis

Conduct eigenanalysis of S (or P; the results will be identical as long as S == P).

[Link] <- eigen(S)

Step 4: Consider the Eigenvalues

There will be as many eigenvalues as there are variables.
[Link]$values %>% round(3) # Eigenvalues
[1] 4.439 1.768 1.531 0.734 0.473 0.362 0.253 0.244 0.143 0.053

The eigenvalues sum to the total of the diagonal of the matrix analyzed. If the data have been
normalized (S) or are based on a correlation matrix (P), the diagonal values are all 1 (see above)
and hence the eigenvalues sum to the number of variables:
sum([Link]$values)
[1] 10

Eigenvalues are generally expressed as a proportion of the total amount of variation:

[Link] <- [Link]$values / sum([Link]$values)
[Link] %>% round(3) # rounding for display
[1] 0.444 0.177 0.153 0.073 0.047 0.036 0.025 0.024 0.014 0.005

The first eigenvalue always accounts for the most variance, and the proportion of variation
declines for each subsequent eigenvalue. The cumulative variance explained is analogous to the
R2 value from a regression. Each eigenvalue is associated with a principal component.

One of the important decisions to be made here is how many eigenvalues to retain; see ‘How
Many Principal Components?’ section below for details. In this example, we will focus on the
first three principal components. Together, these three PCs account for 0.444 + 0.177 + 0.153 =
0.774, or 77% of the variance in the 10 response variables.

Step 5: Consider the Eigenvectors

Every eigenvalue has an associated eigenvector:
[Link]$vectors %>% round(2) # 1 vector per eigenvalue
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] -0.30 -0.49 -0.03 0.12 0.29 0.53 0.06 0.19 0.03 -0.50
[2,] -0.39 -0.20 0.18 0.00 -0.14 -0.42 0.24 0.69 0.13 0.19
[3,] 0.02 -0.35 -0.63 0.07 0.28 -0.54 0.20 -0.24 -0.01 -0.09
[4,] 0.20 -0.29 -0.53 -0.16 -0.63 0.28 -0.16 0.22 -0.03 0.16
[5,] -0.40 0.26 -0.11 0.07 -0.30 0.03 0.37 -0.12 -0.70 -0.17
 [6,] -0.37 0.28 -0.24 0.17 -0.14 0.24 0.36 -0.23 0.64 0.16
[7,] -0.27 0.35 -0.34 0.44 0.17 -0.03 -0.61 0.28 -0.07 0.02
[8,] -0.40 -0.13 0.16 -0.23 -0.39 -0.29 -0.44 -0.33 0.20 -0.40
[9,] -0.37 -0.40 0.12 0.02 0.13 0.13 -0.19 -0.33 -0.21 0.68
[10,] -0.23 0.25 -0.24 -0.82 0.34 0.10 -0.05 0.13 -0.02 0.05

In this matrix of eigenvectors, the rows correspond to the original variables, in order, and the
columns correspond to the eigenvalues. The elements of an eigenvector are known as
the loadings.

The eigenvector associated with each eigenvalue is the coefficients of the linear equation that
‘blends’ the variables together for that axis. Interpretation is based on these loadings. Since the
first three PCs account for most of the variation in this dataset, let’s focus on them. We’ll add
row and column names to aid in interpretation.
loadings <- round([Link]$vectors[,1:3], digits = 3)
rownames(loadings) <- colnames([Link])
colnames(loadings) <- c("PC1", "PC2", "PC3")
loadings
PC1 PC2 PC3
height -0.297 -0.494 -0.026
[Link] -0.390 -0.199 0.178
[Link] 0.016 -0.346 -0.635
[Link] 0.197 -0.290 -0.527
[Link] -0.399 0.257 -0.114
[Link] -0.371 0.284 -0.243
wingsprea -0.272 0.354 -0.344
hoodmass.g -0.396 -0.133 0.158
tubemass.g -0.371 -0.400 0.124
wingmass.g -0.234 0.249 -0.238

When summarizing a PCA in tabular form, two additional elements are often added:

 Proportion of variance explained by each PC

 Cumulative proportion of variance explained by this and larger PCs

An example of this table is below. The proportions of variance explained are helpful to remind
us of the relative importance of the PCs, though these proportions obviously should not be
interpreted like loadings!

PC1 PC2 PC3

 - - -
height 0.29 0.49 0.02
7 4 6

- -
0.17
[Link] 0.39 0.19
8
0 9

- -
0.01
[Link] 0.34 0.63
6
6 5

- -
0.19
[Link] 0.29 0.52
7
0 7

- -
0.25
[Link] 0.39 0.11
7
9 4

- -
0.28
[Link] 0.37 0.24
4
1 3

- -
0.35
wingsprea 0.27 0.34
4
2 4

- -
0.15
hoodmass.g 0.39 0.13
8
6 3
 - -
0.12
tubemass.g 0.37 0.40
4
1 0

- -
0.24
wingmass.g 0.23 0.23
9
4 8

0.44 0.17 0.15

Proportion of variance explained
4 7 3

Cumulative proportion of 0.44 0.62 0.77

variance explained 4 1 4

Interpretation of the loadings is subjective, but is based on the magnitude of the loadings and
their relative direction (which variables have loadings in the same or opposite directions).
Summerville et al. (2006) restricted their attention to loadings > 0.4.

After examining the magnitude and direction of the loadings, here is how we might interpret the
results of this PCA:

 PC1 gives roughly equal weight to all variables except [Link] and [Link], and can
be interpreted as a measure of pitcher size. This PC accounts for 44% of the variance in
the response variables.
 PC2 gives opposing weights to pitcher height and diameter than to the dimensions of the
fishtail appendage ([Link], [Link], wingsprea), and thus can be interpreted
as a measure of pitcher shape. This PC accounts for 18% of the variance in the response
variables.
 PC3 gives most weight to tube and keel diameters, and might be interpreted as a measure
of leaf width. This PC accounts for 15% of the variance in the response variables.

Step 6: Consider Scores of Sample Units

The score or location of each sample unit on each axis is determined by matrix-multiplying the
original data matrix by the eigenvector for that axis.
For example, the scores for all sample units for the first principal component are:
PC1 <- [Link] %*% [Link]$vectors[,1]
PC1 %>% head()
[,1]
[1,] -1.4648658
[2,] 3.2696887
[3,] 1.0026952
[4,] -0.7149071
[5,] 0.9887680
[6,] -1.2136554

The PC scores for all 87 plants for all PCs are calculated in the same way:
[Link] <- [Link] %*% [Link]$vectors
I haven’t displayed [Link] here as it is large (87 rows x 10 columns … right?).

The PC scores are uncorrelated with (orthogonal to) each other:

cor([Link]) %>% round(3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0
[9,] 0 0 0 0 0 0 0 0 1 0
[10,] 0 0 0 0 0 0 0 0 0 1

This is in strong contrast to the correlations among the original variables – see the
objects S and P that we created in step 2 above. This means that we can analyze or use each PC
independently of the others.

The PC scores can also be graphed (see ‘The Biplot: Visualizing a PCA’ section below).
How Many Principal Components?
By definition, there are as many principal components as there were variables in the dataset.
Each principal component is an axis in the ordination space, and has an associated eigenvalue
and eigenvector.

Use Them All … (Data Rotation)

If we retain all of the principal components, we have rotated the data by expressing the location
of each sampling unit in terms of the new synthetic variables created for each principal
component, but we have not altered the dimensionality of the data. In this case, the rotated data
are in the same relative positions as in the original data; all that has changed is the axes along
which the data cloud is arrayed. This is done, for example, following a NMDS ordination
(see upcoming notes on that topic) so that the axes of the ordination solution explain a decreasing
amount of variation.

… Or Focus on a Few (Data Reduction)

When our goal is data reduction, we want to focus on a small number of principal components
while accounting for as much variation as possible. Since the eigenvalues are arranged in
descending order, each accounting for as much of the remaining variance as possible, the first
few PCs will often explain most of the variation. Later PCs account for little variation and
therefore can be dropped with minimal consequences for interpretation.

In addition, since the PCs are orthogonal to one another we can consider them individually. This
is helpful when they capture different ‘facets’ of the responses – see the above interpretation of
pitcher size (PC1), pitcher shape (PC2), and leaf width (PC3). We do not have to decide which
of our measured variables best serves as an index of that facet, but instead can incorporate them
all. We can use this to test whether pitcher size relates to some predictor, without size being
influenced by pitcher shape.

How Many PCs to Focus on?

There are several ways to determine how many principal components (axes) are worthy of
further examination (Peres-Neto et al. 2005). These ways may give different results.

Visually, we can examine a scree plot:

barplot([Link])
 This is different from but similar to the scree plot we saw with cluster analysis – we would be
interested in a sharp bend or change in slope in the ordered eigenvalues.

We could also focus on those principal components that explain more variation than expected.
For example, these could be the eigenvalues that are greater than the mean eigenvalue. This can
be superimposed onto the scree plot:
barplot([Link])
abline(h = mean([Link]),
col = "red")

See below for the broken stick, another way of calculating the amount of variation each
eigenvalue is expected to explain.
 Finally, this decision could be made quantitatively: to focus on as many eigenvalues as necessary
to account for 75% (Borcard et al. 2018) or 90% (Crawley 2007) of the total variation.

In our case, it is reasonable to focus on the first three components. There appears to be a break
in the scree plot below this, each of these eigenvalues explains more variation than expected, and
together they explain 77% of the variance.

PCA Functions in R
While it is helpful to see the steps involved in a PCA, it would be slightly laborious to work
through them each time you did a PCA. R, of course, contains several PCA functions. They
differ slightly in computational method and in the format and contents of the output. The two
that we will review here are stats::princomp() and stats::prcomp(). Others
include vegan::rda() and ade4::[Link]().
stats::princomp()
The princomp() function closely parallels what we did above, using the eigen() function to
calculate the eigenvalues and eigenvectors. Its usage is:
princomp(x,
cor = FALSE,
scores = TRUE,
covmat = NULL,
subset = rep_len(TRUE, nrow([Link](x))),
fix_sign = TRUE,
...
)

The key arguments are:

 x – a matrix or data frame containing the data

 cor – whether to use the correlation matrix (TRUE) or the covariance matrix (FALSE).
Default is the latter. However, we noted above that the variables need to be normalized if
using the covariance matrix. This step is not built into princomp().
 scores – whether to return the score for each observation. Default is to do so (TRUE).
 fix_sign – whether to set the first loading of each PC to be positive (TRUE). This
ensures that analyses are more likely to be comparable among runs. If FALSE, it is
possible for the loadings to have opposing signs from one run to another (note that this
does not change the correlations among the variables but simply whether large scores are
associated with large or small values of each variable (positive and negative signs,
respectively).
This function returns a variety of components, each of which can be indexed and called:

 sdev – the square roots of the eigenvalues of the variance/covariance matrix. Compare
the square of these values to [Link]$values.
 loadings – the matrix of loadings (i.e., eigenvectors). This matrix is easier to view than
the other loading matrices we’ve seen, but it has a few quirks. For example, the default is
to automatically omit elements < 0.1 (a different cutoff can be specified in
the summary() function below). Compare to [Link]$vectors.
 center –the mean that was subtracted from each element during normalization.
 scale – the standard deviation that each element in each column was divided by during
normalization.
 [Link] – the number of observations (sample units)
 scores – the scores for each observation on each principal component. Compare
to [Link] (or, equivalently, scaled data %*% loadings).
 call – the function that was run

Applying this to our example dataset:

[Link] <- princomp([Link], cor = TRUE)

We can use summary() to obtain a summary of the object we created. Behind the scenes, this
function recognizes the class of the object (princomp) and applies pre-defined actions based on
this class. Knowing this information allows us to find help for it, in the general
form [Link](). For example, ?[Link] provides arguments to display
loadings or not (loadings = FALSE) and, if they are displayed, to restrict attention to those
above a specified value (default is cutoff = 0.1). Once we know these arguments, we can of
course adjust them:
summary([Link], loadings = TRUE, cutoff = 0)
Importance of components:
Comp.1 Comp.2 Comp.3 Comp.4
Comp.5
Standard deviation 2.1069368 1.3296011 1.2373322 0.85674840
0.68802564
Proportion of Variance 0.4439183 0.1767839 0.1530991 0.07340178
0.04733793
Cumulative Proportion 0.4439183 0.6207022 0.7738013 0.84720307
0.89454100
Comp.6 Comp.7 Comp.8
Comp.9
Standard deviation 0.60135218 0.50312405 0.49361594
0.37793991
Proportion of Variance 0.03616244 0.02531338 0.02436567
0.01428386
Cumulative Proportion 0.93070344 0.95601682 0.98038249
0.99466635
Comp.10
Standard deviation 0.230946945
Proportion of Variance 0.005333649
Cumulative Proportion 1.000000000

Loadings:
Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
Comp.8
height 0.297 0.494 0.026 0.122 0.285 0.534 0.062
0.185
[Link] 0.390 0.199 -0.178 -0.004 -0.139 -0.418 0.244
0.688
[Link] -0.016 0.346 0.635 0.071 0.282 -0.537 0.196 -
0.238
[Link] -0.197 0.290 0.527 -0.165 -0.626 0.280 -0.157
0.225
[Link] 0.399 -0.257 0.114 0.072 -0.298 0.034 0.372 -
0.121
[Link] 0.371 -0.284 0.243 0.167 -0.141 0.241 0.357 -
0.234
wingsprea 0.272 -0.354 0.344 0.437 0.173 -0.033 -0.614
0.279
hoodmass.g 0.396 0.133 -0.158 -0.232 -0.391 -0.293 -0.439 -
0.334
tubemass.g 0.371 0.400 -0.124 0.024 0.132 0.133 -0.189 -
0.330
wingmass.g 0.234 -0.249 0.238 -0.821 0.345 0.100 -0.047
0.134
 Comp.9 Comp.10
height 0.028 0.497
[Link] 0.126 -0.186
[Link] -0.005 0.094
[Link] -0.030 -0.162
[Link] -0.697 0.166
[Link] 0.638 -0.165
wingsprea -0.070 -0.024
hoodmass.g 0.201 0.404
tubemass.g -0.209 -0.680
wingmass.g -0.022 -0.050

The results shown here should look very similar – they’re the same as we obtained in our step-
by-step approach above.

To obtain a screeplot of the eigenvalues / principal components we can

apply plot() or screeplot() to an object of class princomp:
plot([Link])

The vegan package includes a similarly named screeplot() function that can be applied to
objects of class ‘princomp’ or ‘prcomp’. It includes a ‘broken stick’ argument (bstick) that
compares the variation explained by each PC to a set of ‘ordered random proportions’ – the idea
is that you would retain as many PCs as explain more variation than expected.
stats::prcomp()
 The prcomp() function uses a different approach, singular value decomposition
(the svd() function; see Matrix Algebra notes), to obtain the eigenvalues and eigenvectors.
Manly & Navarro Alberto (2017) note that this algorithm is considered to be more accurate than
that of princomp(), though I have not encountered differences between them. Its usage is:
prcomp(x,
retx = TRUE,
center = TRUE,
scale. = FALSE,
tol = NULL,
rank. = NULL,
...
)

The arguments are:

 x – the matrix or data frame containing the data

 retx – whether to return the rotated variables. Default is to do so (TRUE).
 center – whether to center each variable on zero (i.e., subtract the mean value). Default
is to do so (TRUE).
 scale. – whether to scale each variable to unit variance. Default is to not do so (FALSE),
though scaling is generally recommended as noted in our worked example (and stated in
the help file). Note the ‘.’ in the argument name!
 tol – tolerance; components will be deleted if their standard deviation is less than this
value times the standard deviation of the first component. Default is to include all
variables.
 rank. – number indicating how many principal components to report. Default is to
report them all. All principal components are still calculated; this just limits attention to
those up to the rank specified.

This function returns many of the same components as from princomp(), but calls them by
different names:

 sdev – the standard deviations of the principal components, also known as the square
roots of the eigenvalues of the variance/covariance matrix. Compare squared values
to [Link]$values.
 rotation – the matrix of loadings (i.e., eigenvectors). Compare to [Link]$vectors.
 x – the scores for each observation on each principal component. Compare
to [Link] (or, equivalently, scaled data %*% loadings).
 center – the mean that was subtracted from each element during normalization.
Compare to apply([Link], 2, mean).
 scale – the standard deviation that each element in each column was divided by during
normalization. Compare to apply([Link], 2, sd).

Applying this to our example dataset:

darl.PCA2 <- prcomp([Link], scale. = TRUE) # Note scale. argument

The summary(), print(), plot(), and screeplot() functions behave very similarly when
applied to an object of class prcomp as they do when applied to an object of class princomp as
shown above. For example, the print() function returns the variance explained by ech
principal component along with their loadings:
print(darl.PCA2)
Standard deviations (1, .., p=10):
[1] 2.1069368 1.3296011 1.2373322 0.8567484 0.6880256 0.6013522
0.5031240
[8] 0.4936159 0.3779399 0.2309469

Rotation (n x k) = (10 x 10):

PC1 PC2 PC3 PC4
PC5
height -0.29725091 -0.4941374 0.02594111 -0.122435477
0.2852089
[Link] -0.38993190 -0.1987772 -0.17804342 0.003956204 -
0.1387466
[Link] 0.01633263 -0.3464781 0.63481575 -0.071324686
0.2822023
[Link] 0.19728430 -0.2902694 0.52651716 0.164670782 -
0.6259843
[Link] -0.39911683 0.2574223 0.11432171 -0.071529207 -
0.2984514
[Link] -0.37099329 0.2836524 0.24337230 -0.166545660 -
0.1412989
wingsprea -0.27156026 0.3541451 0.34374791 -0.437068261
0.1732058
hoodmass.g -0.39632173 -0.1326148 -0.15831785 0.231712503 -
0.3906018
tubemass.g -0.37122939 -0.4001026 -0.12442151 -0.024231381
0.1319643
wingmass.g -0.23420022 0.2493993 0.23750541 0.821358528
0.3448007
 PC6 PC7 PC8 PC9
PC10
height -0.53358014 -0.06213397 0.1851593 -0.028165685
0.49680679
[Link] 0.41785432 -0.24376093 0.6877145 -0.125773703 -
0.18623421
[Link] 0.53670782 -0.19570472 -0.2384156 0.005035915
0.09352417
[Link] -0.28028754 0.15678491 0.2246203 0.030018781 -
0.16165211
[Link] -0.03412138 -0.37171835 -0.1206353 0.696871529
0.16630928
[Link] -0.24090857 -0.35700239 -0.2344052 -0.638290495 -
0.16469409
wingsprea 0.03256432 0.61420874 0.2791219 0.069750311 -
0.02361776
hoodmass.g 0.29328788 0.43946368 -0.3340370 -0.200905130
0.40364300
tubemass.g -0.13345583 0.18870304 -0.3296105 0.209051028 -
0.68033777
wingmass.g -0.09991858 0.04715240 0.1335554 0.022049695 -
0.04993778

Using PC Scores
We can use principal components in several ways during subsequent analyses. When doing so,
however, we must explicitly state which variables were included in the PCA and whether we
used all of the loadings (and, if not, how we decided to include or exclude them). We must also
keep in mind that the interpretation of the results may depend on how we are interpreting those
loadings.

Reduce Dimensionality of Explanatory Variables

One way to use PCs is to apply a PCA to highly correlated explanatory variables and
thereby summarize many explanatory variables in a much smaller number of uncorrelated
 variables (ideally, one or two) that can then be used as explanatory variables in subsequent
analyses. Haugo et al. (2011) provide an example of this approach. One potential criticism of
this approach is that the PC is more abstract than the set of variables from which it is calculated.

Analyze Each PC Independently

Another way to use PCs is to apply a PCA to highly correlated response variables and
then analyze the response as represented by individual PCs. Because the PCs are
uncorrelated with one another, it is appropriate to analyze them in separate univariate hypothesis
tests.

For example, based on our examination of the loadings we interpreted PC1 to be a measure of
pitcher size. We can use PERMANOVA to test for a difference in pitcher size among sites. To
do so, we’ll first combine the scores with the site identity of each plant:
[Link] <- [Link]([Link]$scores,
site = darl$site)
adonis2([Link]$Comp.1 ~ site,
data = [Link],
method = "euc")
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = [Link]$Comp.1 ~ site, data =

[Link], method = "euc")
Df SumOfSqs R2 F Pr(>F)
site 3 100.53 0.26031 9.7362 0.001 ***
Residual 83 285.68 0.73969
Total 86 386.21 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This test indicates that pitcher size differs among the sites. We could follow this test with
pairwise contrasts to determine exactly which sites differ in ‘pitcher size’.

We could also see these differences visually:

ggplot(data = [Link], aes(x = site, y = Comp.1)) +
geom_boxplot() +
geom_jitter(aes(colour = site), width = 0.3, height = 0) +
theme_bw()

Since PCs are orthogonal, we could conduct identical tests for differences among sites with
respect to PC2 (pitcher shape) and PC3 (leaf width). Tests of one PC will always be
independent of tests of other PCs produced by that PCA.

The Biplot: Visualizing a PCA

Base R Graphing
PCs are often examined visually. The first two principal components are often graphed in
a biplot:
biplot([Link])
This kind of graphic is quite common in ordinations. Some comments about it:

 We focus on the first two PCs because, by definition, they explain more of the variation
than any other PCs. To create this in our step-by-step introduction to PCA, we simply
needed to graph the first two PCs: plot([Link])
 Each individual data point (row in the sample × species matrix) is represented by its row
number. We could customize this plot to display them in other ways, such as by having
different symbol shapes and colors for different levels of a grouping variable.
 The scores on each PC axis are centered on zero because of the normalizing that occurred
during the PCA.
 The axes on the bottom and left-hand side of the graph show the values of the first two
PCs. It is common practice to report the amount of variation explained by each axis in
the axis title. For example, this code would change the axis labels of the
biplot: biplot([Link], xlab= "PC1 (44.4%)", ylab= "PC2 (17.7%)")
 A second coordinate system is shown on the top and right-hand side. According to
Everitt & Hothorn (2006, p. 222), this displays the first two loadings associated with
these variables – though I haven’t verified this.
 Each response variable is displayed as a vector (arrow) showing the direction in which it
is most strongly linearly correlated. These correlations are reflected in the loadings for
that variable with each PC. In particular:
o Variables that are perfectly correlated with an axis are parallel to it
o Variables that are uncorrelated with an axis are perpendicular to it
  The angle between vectors shows the strength and direction of the correlation between
those variables. You should be able to find examples of each of these cases in the above
biplot:
o Vectors that are very close to one another represent variables that are strongly
positively correlated
o Vectors that are perpendicular represent variables that are uncorrelated
o Vectors that point in opposite directions represent variables that are strongly
negatively correlated

The ggbiplot Package

The ggplot2 package does not have a built-in capability to draw biplots, but
the ggbiplot package allows their creation using the ggplot2 approach. It is available through
github:
devtools::install_github("vqv/ggbiplot")
library(ggbiplot)
ggbiplot([Link])

Compare to the above biplot produced with base R plotting capabilities.

This function can be combined with other standard ggplot2 functions. For example:
ggbiplot([Link]) +
geom_point(aes(colour = darl$site)) +
labs(colour = "Site", x = "PC1(Size; 44.4%)", y = "PC2 (Shape; 17.7%)") +
theme_bw()
I’ve changed the theme and the axis labels, and color-coded each plant by the site it came from.
However, remember that site identity was not part of the PCA.

Conclusions
Since PCA identifies the best linear fit, it therefore is, logically, intended for data with linear
relationships among variables. PCA is appropriate under certain scenarios:

 if the sample units span a short gradient (i.e., beta diversity is low; little turnover among
sample units)
 if the data matrix is not sparse (i.e., does not contain many nonzero values)
 when variables (attributes) are highly correlated with one another

Examples where PCA is often appropriate include:

 LiDAR data (multiple highly correlated and continuously distributed variables)

 fuels data
 morphological data
 physiological data

PCA has been used for community-level data (i.e., sample unit × species matrices) in the past,
but these applications are considered unsuitable (Beals 1971; Minchin 1987; Legendre &
Legendre 2012). Problems include:

 Assumes continuously distributed data (no nominal or ordinal categories)

  Assumes multivariate normality. This is particularly important if you are going to use a
PCA to make statistical inferences (vs. describing a data set).
 PCA works best with highly correlated response variables, but the abundances of many
species are only weakly correlated with one another.
 Sensitive to outliers. If data are highly skewed, the first few principal components will
separate the extreme values from the rest of the sample units.
 “The quality of ordination by PCA is completely dependent on how well relationships
among the variables can be represented by straight lines” (McCune & Grace 2002, p.
116).
 Need more sample units than variables measured (technicality; I’m not sure how
important this really is).
 Implicitly uses the Euclidean distance measure
 May yield horseshoe effect (artificial curvature) in second and higher axes if sample units
span a long gradient. This occurs because shared zeroes are interpreted as an indicator of
a positive relationship, and because two species may be positively related to each other at
some points and negatively related at other points. See Fig. 19.2 in McCune & Grace
(2002) for an example of this.
 ‘Axis-centric’? (term borrowed from Mike Kearsley at Northern Arizona University) –
we tend to assign more meaning to axes than may be warranted.
 Tendency to interpret proximity of apices of variables (vectors) rather than the angles
separating them.

Although PCA is inappropriate for analyzing community-level data, there are certain
applications in which it is appropriate even with these data. In particular, as we’ll soon see, PCA
is often used after a NMDS ordination to express the ordination solution in such a way that as
much variation as possible is explained by the first axes. When doing so, all principal
components are retained and thus the dimensionality of the ordination solution is unchanged; the
rotated data are in the same relative positions as in the original ordination solution.

References
Beals, E.W. 1971. Ordination: mathematical elegance and ecological naivete. Journal of
Ecology 61:23-35.

Borcard, D., F. Gillet, and P. Legendre. 2018. Numerical ecology with R. 2nd edition. Springer,
New York, NY.

Crawley, M.J. 2007. The R book. John Wiley & Sons, Hoboken, NJ.

Everitt, B.S., and T. Hothorn. 2006. A handbook of statistical analyses using R. Chapman &
Hall/CRC, Boca Raton, LA.

Goodall, D.W. 1954. Objective methods for the classification of vegetation. III. An essay in the
use of factor analysis. Australian Journal of Botany 2:304-324.
Gotelli, N.J., and A.M. Ellison. 2004. A primer of ecological statistics. Sinauer, Sunderland,
MA.

Haugo, R.D., C.B. Halpern, and J.D. Bakker. 2011. Landscape context and tree influences shape
the long-term dynamics of forest-meadow ecotones in the central Cascade Range,
Oregon. Ecosphere 2:art91.

Hotelling, H. 1933. Analysis of a complex of statistical variables in principal

components. Journal of Experimental Psychology 24:417-441, 493-520.

Legendre, P., and L. Legendre. 2012. Numerical ecology. 3rd English edition. Elsevier,
Amsterdam, The Netherlands.

Manly, B.F.J., and J.A. Navarro Alberto. 2017. Multivariate statistical methods: a primer.
Fourth edition. CRC Press, Boca Raton, FL.

McCune, B., and J.B. Grace. 2002. Analysis of ecological communities. MjM Software Design,
Gleneden Beach, OR.

Minchin, P.R. 1987. An evaluation of the relative robustness of techniques for ecological
ordination. Vegetatio 69:89-107.

Pearson, K. 1901. On lines and planes of closest fit to systems of points in space. Philosophical
Magazine, Sixth Series 2:559-572.

Peres-Neto, P.R., D.A. Jackson, and K.M. Somers. 2005. How many principal components?
stopping rules for determining the number of non-trivial axes revisited. Computational Statistics
& Data Analysis 49:974-997.

Summerville, K.S., C.J. Conoan, and R.M. Steichen. 2006. Species traits as predictors of
lepidopteran composition in restored and remnant tallgrass prairies. Ecological
Applications 16:891-900.

Media Attributions
 Darlingtonia.californica_Jepson
 PCA.barplot1
 PCA.barplot2
 [Link]
 [Link]
 [Link]
 [Link]
 [Link].gg2

Previous/next navigation
Previous: Overview of Ordination Methods
 Next: NMDS
BACK TO TOP
LICENSE

Applied Multivariate Statistics in R Copyright © 2023 by Jonathan D.

Bakker is licensed under a Creative Commons Attribution-NonCommercial-
ShareAlike 4.0 International License, except where otherwise noted.

SHARE THIS BOOK

Powered by Pressbooks

 Guides and Tutorials

 |Pressbooks Directory
 |Contact

Pressbooks on YouTubePressbooks on Twitter