Biplot using base graphic functions in R

See Video ⮞ ☝

AGRON Stats

June 26, 2018

• Introduction
• Import data set
• Principal Component analysis
• Display biplot
• Add axis, title and labels
• Add points
• Add ellipse
• Add axis for variable vectors
• Add labels to second axis
• Add box and ablines
• Add variable vectors
• Add lables to variable vectors
• Add unit variance circle
• Adding Legend

Introduction

This tutorial builds a principal component biplot with base R graphics using measurements from genuine and counterfeit banknotes. A biplot displays sample scores and variable information for selected components, allowing an exploratory view of structure in the data.

The example uses the banknote data introduced in earlier PCA tutorials. Confirm the measurement definitions and group labels in the source file before interpreting the plotted clusters.

Also go through some video tutorials to understand the data set, principal component analysis and biplot interpretation — PCA_R & Biplot_PCA_R.

We will construct and customize the score display with base graphics, then overlay variable vectors on a second scale. The result is flexible, but the two plotted scales must be interpreted separately.

The workflow uses base R plotting functions to add the score points, group ellipses, axes, labels, and variable vectors in layers.

Clear the workspace only if necessary and save any objects you need to retain. Closing graphics devices and clearing the console are optional; neither step changes the data.

rm(list = ls(all = TRUE))
graphics.off()
shell("cls")

Import data set

Import the CSV file using the correct path. A project-relative path improves reproducibility when the analysis is moved to another computer.

Set header = TRUE when the first row contains column names, then inspect the imported data with head() and str(). Confirm that measurements are numeric and the group identifier is categorical.

data <- read.csv(file = "biplot.csv", 
                 header = TRUE)
head(data)

#      type Length  Left Right Bottom  Top Diagonal
# 1 genuine  214.4 130.1 130.3    9.7 11.7    139.8
# 2 genuine  214.9 130.5 130.2   11.0 11.5    139.5
# 3 genuine  214.9 130.3 130.1    8.7 11.7    140.2
# 4 genuine  215.0 130.4 130.6    9.9 10.9    140.3
# 5 genuine  214.7 130.2 130.3   11.8 10.9    139.7
# 6 genuine  215.0 130.2 130.2   10.6 10.7    139.9

Principal Component analysis

Fit PCA to the numeric measurements. In base R, prcomp() centers columns by default; specify whether to scale them according to the measurement scales and scientific question.

The input to PCA should contain the features of interest, excluding identifiers and group labels. The expression data[, -1] removes the first column only; verify that this is the intended grouping column before using it.

Centering subtracts each variable mean. Scaling divides centered columns by their standard deviations, giving each variable unit variance. Use scaling when the variables’ units or ranges would otherwise dominate the analysis; same units alone do not always justify omitting it.

Inspect the fitted object and its summary for component standard deviations, loadings, and proportions of variance. Check that the selected components are interpretable for the study purpose.

require(stats)
pc <- prcomp(x = data[ , -1],
             center = TRUE, 
             scale. = FALSE)
print(pc)
summary(pc)

# Standard deviations (1, .., p=6):
# [1] 1.7321388 0.9672748 0.4933697 0.4412015 0.2919107 0.1884534
# 
# Rotation (n x k) = (6 x 6):
#                  PC1         PC2        PC3        PC4         PC5         PC6
# Length    0.04377427  0.01070966 -0.3263165 -0.5616918  0.75257278 -0.09809807
# Left     -0.11216159  0.07144697 -0.2589614 -0.4554588 -0.34680082  0.76651197
# Right    -0.13919062  0.06628208 -0.3447327 -0.4153296 -0.53465173 -0.63169678
# Bottom   -0.76830499 -0.56307225 -0.2180222  0.1861082  0.09996771  0.02221711
# Top      -0.20176610  0.65928988 -0.5566857  0.4506985  0.10190229  0.03485874
# Diagonal  0.57890193 -0.48854255 -0.5917628  0.2584483 -0.08445895  0.04567946
# Importance of components:
#                           PC1    PC2     PC3     PC4     PC5    PC6
# Standard deviation     1.7321 0.9673 0.49337 0.44120 0.29191 0.1885
# Proportion of Variance 0.6675 0.2082 0.05416 0.04331 0.01896 0.0079
# Cumulative Proportion  0.6675 0.8757 0.92983 0.97314 0.99210 1.0000

Display biplot

The tutorial focuses on the first two components. Calculate their cumulative variance from the fitted object rather than relying on a hard-coded percentage; the value depends on the preprocessing and feature set.

Base graphics are available in a standard R installation. The following commands set up the plotting device and add each graphical layer.

Use par() to inspect current graphics parameters or set them for subsequent plotting calls. Record any non-default settings needed to reproduce the figure.

require(graphics)
par()

Set shared parameters once with par() when several low-level plotting functions will use the same style. Restore previous settings after plotting if the change should not affect later figures.

Choose a square plotting region with pty = "s" when equal axis scaling is important. Adjust text size and font only as needed for legibility; avoid shrinking labels to fit an overcrowded plot.

Set title and axis-label colors consistently. Ensure that contrast remains strong when the figure is printed or viewed by readers with color-vision differences.

The foreground color parameter can also control axes and plot borders. Use it consistently so that axes do not compete visually with the data.

The las setting controls axis-label orientation. Horizontal labels can improve readability when category or value labels are long.

require(graphics)
par(pty = "s",
    cex.main = 1.2,
    cex.lab = 1,
    font.main = 2,
    font.lab = 2,
    family = "sans",
    col.main = "gray10",
    col.lab = "gray10",
    fg = "gray10",
    las = 1)

Start a new plotting frame with plot.new() when assembling a custom figure from low-level functions.

Set the plotting limits from the score coordinates for the selected components. Include enough padding for point labels and group ellipses, and retain comparable aspect scaling when interpreting distances.

Use plot.window() to define x- and y-limits and the aspect ratio. Because the next layer uses a different scale for variable vectors, make the transformation between score and loading coordinates explicit.

plot.new()
plot.window(xlim = c(-4, 4), 
            ylim = c(-4, 4), 
            asp = 1)

Add axis, title and labels

Use axis() to draw the required axes, specifying the side, tick positions, and labels deliberately. Check that each axis corresponds to the plotted coordinate scale.

In base graphics, side = 2 places the y-axis on the left. Set tick locations and labels to match the score coordinates and avoid implying precision beyond the plotted scale.

Add a concise title and axis labels with title(). Use alignment and margin settings to keep text readable and consistent with the plotted component scores.

Label the axes with the component names and their variance proportions. Obtain those proportions directly from the PCA summary and state whether the PCA was based on centered or standardized variables.

For PC1, the proportion of variance is in the “Proportion of Variance” row and PC1 column of the importance matrix. Use explicit row and column indexing to avoid selecting an unintended value.

Repeat the same calculation for PC2 and label it clearly. Confirm that the two percentages are consistent with the fitted PCA output.

Place axis labels in the figure margins with suitable line and alignment settings. Preview the result to ensure labels are not clipped by the graphics device.

axis(side = 1, 
     at = c(-4, -2, 0, 2, 4),
     labels = TRUE)
axis(side = 2, 
     at = c(-4, -2, 0, 2, 4),
     labels = TRUE)
title(main = "Biplot for PCs of bank note data", 
      line = 3,
      adj = 0)
title(xlab = paste("PC 1 (", 
                   round(summary(pc)$importance[2]*100, 
                         digits = 1),
                   "%)", 
                   sep = ""), 
      ylab = paste("PC 2 (", 
                   round(summary(pc)$importance[5]*100, 
                         digits = 1),
                   "%)", 
                   sep = ""), 
      line = 2,
      adj = 0.5)

Add observations using their score coordinates for the selected components. Scores are projected coordinates, not the original measurements; keep axis labels explicit about the components displayed.

Add points

Use the group variable to map point symbols or colors. Do not assign group membership from row order unless the source data guarantee that order; derive the mapping from the actual group column.

The example uses two groups of banknotes. Verify the count and labels for each group in the imported data, then map symbols consistently using the factor levels.

points(x = pc$x[,1:2],
       pch = c(rep(16, times = 100),
               rep(17, times = 100)),
       cex = 1,
       col = c(rep("darkcyan", times=100),
               rep("orangered", times=100)))

The ellipse package can draw group ellipses from score coordinates. Install and load it as required, and verify that the ellipse function’s input and assumptions match the data.

Add ellipse

Construct each ellipse from the score coordinates for one group. Ensure that the covariance matrix and center are calculated from the same selected observations and components.

The center is the group mean in score space. The level argument controls the probability contour under the ellipse model; it is not automatically a confidence interval for the group mean. Interpret such contours as distributional summaries only when assumptions are reasonable.

Use polygon settings to distinguish group contours without obscuring the observations. Keep transparency and line weight appropriate for the final display size.

Overlaying variable vectors requires a second coordinate scale or a carefully defined biplot transformation. par(new = TRUE) starts a new plotting layer; it does not by itself make score and loading scales directly comparable.

Set limits for the variable-vector layer separately and document the scaling transformation used to overlay it on the score plot.

library(ellipse)
polygon(ellipse(x = cor(pc$x[1:100, 1], pc$x[1:100, 2]),
                centre = colMeans(pc$x[1:100, 1:2]),
                level = 0.85),
        border = "gray10",
        lty = "solid",
        lwd = 1.5,
        col = adjustcolor("darkcyan", alpha.f = 0.20))
polygon(ellipse(x = cor(pc$x[101:200, 1], pc$x[101:200, 2]),
                centre = colMeans(pc$x[101:200, 1:2]),
                level = 0.85),
        border = "gray10",
        lty = "solid",
        lwd = 1.5,
        col = adjustcolor("orangered", alpha.f = 0.20))

Draw the second set of axes with ticks and labels that correspond to the variable-vector coordinates. Make the distinct score and loading scales clear to the reader.

Add axis for variable vectors

Add the right-side y-axis only if the variable-vector scale requires it, and label it to distinguish it from the score axis.

Use mtext() to add labels in plot margins. Choose the side and line position carefully to avoid overlap with the plot or clipping.

Add the second vertical-axis label only when it clarifies the loading scale. Keep wording distinct from the score-axis label.

The las parameter controls the orientation of axis text. Choose an orientation that remains readable without crowding the plotting region.

par(new = TRUE, las = 1)
plot.window(xlim = c(-1, 1), 
            ylim = c(-1, 1), 
            asp = 1)
axis(side = 3, 
     at = c(-1, 0.5, 0, -0.5, 1), 
     labels = TRUE, 
     col = "navy",
     col.ticks = NULL,
     lwd = 2,
     col.axis = "navy")
axis(side = 4, 
     at = c(-1, 0.5, 0, -0.5, 1), 
     labels = TRUE, 
     col = "navy",
     col.ticks = NULL,
     lwd = 2,
     col.axis = "navy")

Add a border and zero reference lines with box() and abline() where useful. Reference lines should clarify the origin without overpowering the score pattern or variable vectors.

Add labels to second axis

If you want to add labels to second axis then use mtext() function. It will write text into the margins of a plot. In text argument you can specify the character or expression vector specifying the text to be written. The value three in side argument specify the placement of text will take place in top margins of the figure region. Again you can specify the cex, font, font family, col and pos (line position) for the axis label.

Similarly use the same function for adding label to the Y axis for second PC. Four value for side argument side = 4 specify the placement of the axis label will take place on the right side of the margins of the plot.

I have used an additional las argument. The value three las = 3 specify the style of the Y axis label will be vertical.

mtext((text = "PC 1 rotations"), 
      side = 3,
      cex = 1,
      font = 2,
      family = "sans",
      col = "gray10", 
      line = 2) 
mtext((text = "PC 2 rotations"), 
      side = 4,
      cex = 1,
      font = 2,
      family = "sans",
      col = "gray10", 
      line = 2,
      las = 3)

Add box and ablines

Now let’s add box and lines to the plot. By simply using box() function will add box around the plot. To draw horizontal and vertical lines for PC1 and PC2 use abline() function. Set zero value for both vertical v = 0 and horizontal lines h = 0. Setting value two for lty argument will draw dashed lines. You can also specify color for the lines in col argument.

box()
abline(v = 0, h = 0, lty = 2, col = "grey25")