You can code it by hand…

# scale and center table
d <- as.data.frame(scale(iris[, 1:4], center = TRUE, scale = TRUE))

# covariance matrix
covariance_matrix <- cov(d)

# eigenvalues
lambdas <- eigen(covariance_matrix)$values
importance_principal_components <- lambdas / sum(lambdas)

# eigenvectors
v <- eigen(covariance_matrix)$vectors
principal_components <- t(t(v) / colSums(v))

… but you don’t need to in R (see next slide).

You can see information about the principal components (PCs) with the summary() function.

summary(pca)

Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

Ignore the row Standard deviation for now.

• The row Proportion of Variance is the variance explained by each PC.
• The row Cumulative Proportion is the proportion of variance explained by the PCs + the proportion explained by PCs with lower index. For instance the cumulative proportion of PC3 is proportion of PC1 + proportion of PC1 + proportion of PC3.
• You want to select PCs that, together, explain circa 70-80 % of the variance, at least.
• This means that you want the PCs for which Cumulative Proportion \(\geq 0.70\) or \(\geq 0.80\).
• In this example, I am happy selecting only the first two PCs (Cumulative Proportion = 0.96).

You can also see the results before by plotting columns against each other. Usually, PCA is to avoid doing this, but here we will do it for understanding PCA.

You can see that our interpretation from PCA is correct.

Principal Coordinate Analysis (PCoA), also known as metric Multi-Dimensional Scaling (mMDS) is similar to PCA, but it needs distance matrices.

For example, the Euclidean distance: \(D = \sqrt{\sum{(x_i - x_j) ^ 2}}\)

PCoA tries to order data on a plot, where the distance between points is proportional to the distance between data.

In this figure, if difference in color indicates the dissimilarity between points, then we want points with more similar color closer together.

You can also see the distance matrix using image(). Unfortunately, image() rotates the matrix 90 degrees, i.e. the rows and column of the image are inverted (but it does not matter, as the matrix is symmetric).

I use hcl.colors() to change the color palette.

In base R, cmdscale(d = <distance>) performs a PCoA, where <distance> is the distance matrix (calculated before).

distance <- dist(iris[, 1:4])
pcoa <- cmdscale(distance)

plot(pcoa)

Are PCoA results much different from PCA?

par(mfrow = c(1, 2), mar = c(4, 4, 2, 2))
pcoa <- cmdscale(distance, k = 2)
plot(pcoa, col = color_for_species(iris$Species))
plot(pca$x[, 1], -pca$x[, 2], col = color_for_species(iris$Species))

Not in this example, but in general they can be.

There are several nMDS algorithms in R, but none in base R. I will use isoMDS() from the MASS package, simply because MASS is often installed together with other packages.

I need to add a small amount to the distance matrix (1e-9 = 0.000000001) because in MASS::isoMDS() distances between rows cannot be zero. Two rows in iris are identical, therefore they have distance = 0. In general, you don’t need to do this.

nMDS <- MASS::isoMDS(distance + 1e-9)  # distances cannot be zero

initial  value 3.025865 
iter   5 value 2.637651
final  value 2.582478 
converged

The algorithm works through iteration and stops when a value (called “stress”) reaches a certain relative value. Don’t think too much about it, just remember that nMDS moves points around several times until it decides that are good enough. In this example, isoMDS() moves around points 5 tiems (5 iterations) before stopping.

Are nMDS results much different from PCoA?

par(mfrow = c(1, 2), mar = c(4, 4, 2, 2))
plot(pcoa, col = color_for_species(iris$Species))
plot(nMDS$points, col = color_for_species(iris$Species))

Not in this example, but in general they can be.

• Partition \(n\) observations into \(k\) clusters (groups).
• Often done after PCA or PCoA.
• Useful to find groups in your data.
• In base R, use kmeans().
• Out of the scope of this course, but useful to remember: if you want to identify groups in your data, do a PCA (prcomp()) or PCoA (cmdscale()) and then a K-mean clustering (kmeans()).

In the figure above, colors show the clusters and shapes the species.