It all started with a comment to always scale the input variables before doing principal components analysis....
The question asks why the PCA biplots generated with stats::biplot.prcomp (in base R) and factoextra::fviz_pca_biplot (built on ggplot2) "look different". It turns out that the plots differ in two ways:
- biplot.prcomp plots the principal components while fviz_pca_biplot plots the principal components scaled to have unit variance.
- biplot.prcomp and fviz_pca_biplot use different x:y aspect ratio to overlay principal components (points) and loadings (arrows) in an aesthetically pleasing way.
Okay, let's do principle component analysis. Don't forget to scale the data beforehand or to set scale = TRUE in the prcomp call. None of what follows holds if the data $\mathbf{X}$ is not standardized.
X <- iris[, 1:4]
X <- scale(X)
n <- nrow(X)
pca <- prcomp(X, scale = TRUE)
Recall that principal component analysis decomposes the matrix $\mathbf{X}$ into three components:
$$
\mathbf{X} = \mathbf{U}\mathbf{S}\mathbf{V}^\top
$$
where $\mathbf{V}$ are the principal axes (the eigenvectors), $\mathbf{S}$ are the standard deviations of the principal axes (the square roots of the eigenvalues) and $\mathbf{U}\mathbf{S}$ are the principal components (aka scores). The standardized principal components $\mathbf{U}$ have unit variance. The loadings are defined as $\mathbf{V}\mathbf{S}/\sqrt{n-1}$.
# the standard deviations of the principal components (the square roots of the eigenvalues)
S <- pca$sdev * sqrt(n - 1) # same as sqrt(eigen(t(X) %*% X, symmetric = TRUE)$values)
the principal axes (the eigenvectors)
V <- pca$rotation # same as eigen(t(X) %*% X, symmetric = TRUE)$vectors
the principal components (the observations projected on the principal axes)
US <- pca$x # same as X %*% V
scores <- US
scores_unit_var <- divide(US, S) * sqrt(n-1)
loadings <- multiply(V, S) / sqrt(n - 1)
@amoeba explains all this in great detail in Positioning the arrows on a PCA biplot
biplot.prcomp plots the "raw" principal components.
# Compute the x:y ratio to overlay PCs and loadings into one plot
# See source code for how biplot.prcomp computes the ratio
# https://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/biplot.R
xy_ratio <- default_xy_ratio(scores, loadings)
biplot_by_hand(scores, loadings, xy_ratio)
fviz_pca_biplot plots the standardized principal components.
# Compute the x:y ratio to overlay standardized PCs and loadings into one plot
# See source code for how factominer computes the ratio
# https://rdrr.io/cran/factoextra/src/R/fviz_pca.R
xy_ratio <- factominer_xy_ratio(scores_unit_var, loadings)
biplot_by_hand(scores_unit_var, loadings, xy_ratio / sqrt(n - 1))
Complete R code listing to reproduce all figures. I use ggplot2.
library("factoextra")
library("tidyverse")
plot_pcs <- function() {
ggplot(
mapping = aes(PC1, PC2)
)
}
add_points <- function(p, data) {
p +
geom_point(
shape = "o",
size = 3,
data = as_tibble(data)
)
}
add_arrows <- function(p, data) {
p +
geom_segment(
aes(
x = 0, xend = PC1,
y = 0, yend = PC2
),
inherit.aes = FALSE,
data = as_tibble(data),
color = "red",
arrow = arrow(length = unit(0.1, "cm"))
)
}
fix_aspect_ratio <- function(p) {
layer <- p$layers[[1]]
xy_limits <- range(layer$data[, unlist(p$labels)])
p +
scale_x_continuous(
limits = xy_limits
) +
scale_y_continuous(
limits = xy_limits
) +
coord_fixed()
}
biplot_by_hand <- function(x, y, xy_ratio) {
plot_pcs() %>%
add_points(x) %>%
add_arrows(y / xy_ratio) %>%
fix_aspect_ratio()
}
multiply <- function(mat, vec) {
x <- mat %*% diag(vec)
dimnames(x) <- dimnames(mat)
x
}
divide <- function(mat, vec) {
multiply(mat, 1 / vec)
}
X <- iris[, 1:4]
X <- scale(X)
n <- nrow(X)
pca <- prcomp(X, scale = TRUE)
the standard deviations of the principal components (the square roots of the eigenvalues)
S <- pca$sdev * sqrt(n - 1) # same as sqrt(eigen(t(X) %*% X, symmetric = TRUE)$values)
the principal axes (the eigenvectors)
V <- pca$rotation # same as eigen(t(X) %*% X, symmetric = TRUE)$vectors
the principal components (the observations projected on the principal axes)
US <- pca$x # same as X %*% V
scores <- US
scores_unit_var <- divide(US, S) * sqrt(n - 1)
loadings <- multiply(V, S) / sqrt(n - 1)
Compute the x:y ratio to overlay PCs and loadings in the same plot
See source code for how biplot.prcomp computes the ratio
https://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/biplot.R
default_xy_ratio <- function(x, y) {
unsigned.range <- function(x) {
c(-abs(min(x, na.rm = TRUE)), abs(max(x, na.rm = TRUE)))
}
rangx1 <- unsigned.range(x[, 1L])
rangx2 <- unsigned.range(x[, 2L])
rangy1 <- unsigned.range(y[, 1L])
rangy2 <- unsigned.range(y[, 2L])
max(rangy1 / rangx1, rangy2 / rangx2)
}
xy_ratio <- default_xy_ratio(scores, loadings)
biplot_by_hand(scores, loadings, xy_ratio)
Compute the x:y ratio to overlay standardized PCs and loadings in the same plot
See source code for how factominer computes the ratio
https://rdrr.io/cran/factoextra/src/R/fviz_pca.R
factominer_xy_ratio <- function(scores, loadings, scale. = 0.5) {
r <- min(
(max(scores[, 1L]) - min(scores[, 1L]) / (max(loadings[, 1L]) - min(loadings[, 1L]))),
(max(scores[, 2L]) - min(scores[, 2L]) / (max(loadings[, 2L]) - min(loadings[, 2L])))
)
r / scale.
}
xy_ratio <- factominer_xy_ratio(scores_unit_var, loadings)
biplot_by_hand(scores_unit_var, loadings, xy_ratio / sqrt(n - 1))
scale = TRUEinprcomp? – dipetkov Jun 29 '22 at 22:20?biplotproduces the advice, "There are many variations on biplots (see the references)..." which at least heads you in a productive direction. Indeed, one of the references is an entire book on this plot (!), "J.C. Gower and D. J. Hand (1996). Biplots. Chapman & Hall." A search for it turns up another more recent book, https://onlinelibrary.wiley.com/doi/book/10.1002/9780470973196. – whuber Jun 30 '22 at 16:48