PCA of iris

The following is a very brief explanation of how the PCA rotates data to maximise its variance along orthogonal axes. It intentionally does not touch on the dimensionality reduction aspect.

First we create a small plotting helper that we can reuse with new data; we’ll use coord_fixed to better illustrate the variability across the two dimensions.

plot = function (data) {
    data = data %>%
        as.data.frame() %>%
        mutate(Species = iris$Species)
    cols = as.list(setNames(colnames(data), c('x', 'y', 'color')))
    ggplot(data, do.call(aes_string, cols)) + geom_point() + coord_fixed()
}

To see the data before its transformation, let’s plot the original iris distribution of petal width and length:

(petal_data = select(iris, Petal.Width, Petal.Length))

plot(petal_data)

Now rotate the data by performing PCA:

pca = prcomp(petal_data)
(rotated_petals = as.data.frame(pca$x))

And plot the rotated data:

plot(rotated_petals)

The prcomp function also gives us the rotation matrix \(r\) of the PCA, which is a matrix of the rotation’s eigenvectors:

pca$rotation

                   PC1        PC2
Petal.Width  0.3877188 -0.9217777
Petal.Length 0.9217777  0.3877188

(Incidentally, this corresponds to a rotation angle \(\theta = {\mathop{\rm acos}r_{11}} \approx 67°\).)

We can apply the rotation ourselves via simple matrix multiplication:

rotated_petals2 = as.matrix(petal_data) %*% pca$rotation
plot(rotated_petals2)

Note that this figure looks almost identical to rotated_petals; the only difference is the \(x\)- and \(y\)-axis labelling: by default, prcomp centers the data before rotation; If we wanted to do the same, we’d have to column-wise subtract pca$center from our rotated data (which is identical to colMeans(petal_data)).

Finally we can visualise the rotations’ eigenvectors and eigenvalues (which are both factors of the loadings) in a so-called biplot. pca$sdev, the rotation’s standard deviation, is the square root of the eigenvalues; we therefore get:

(loadings = sweep(pca$rotation, 2, pca$sdev, `*`))

                   PC1         PC2
Petal.Width  0.7418752 -0.17500689
Petal.Length 1.7637628  0.07361153

A biplot is the PCA plot augmented by arrows corresponding to the loadings:

plot(rotated_petals) + arrows(loadings)

If, instead, we plot the eigenvectors directly, we can see immediately how the data has been rotated from its original form into the PCA:

plot(rotated_petals) + arrows(pca$rotation)

The arrows do indeed correspond to the original (centred) axes of the iris petal data, and the plot’s clockwise tilt by \(67°\) becomes apparent.