Multivariate data analysis involves the simultaneous observation and analysis of more than two outcome variables for each individual or object. Because human vision is naturally limited to two or three dimensions, visualizing high-dimensional data requires clever projection, reduction, and encoding techniques.

This note covers foundation and advanced visualization techniques for multivariate data, complete with reproducible R code and detailed interpretations.

1. Introduction & Dataset Setup

To demonstrate these techniques, we will use the built-in iris dataset, which contains 150 cases of iris flowers characterized by four numeric variables:

Let’s load the necessary libraries and inspect the data structure.

# Load required libraries
library(ggplot2)
library(GGally)
library(reshape2)
library(lattice)

# Inspect the dataset
data(iris)
head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

2. Scatterplot Matrices (Pairs Plots)

A scatterplot matrix displays pairwise relationships between all variables in a dataset. For \(p\) variables, it produces a \(p \times p\) grid of plots where the off-diagonal cells are bivariate scatterplots, and the diagonal cells show univariate distributions (histograms or density plots).

# Using GGally for a publication-ready pairs plot colored by species
ggpairs(iris,
        columns = 1:4,
        aes(color = Species, alpha = 0.5),
        lower = list(continuous = wrap("smooth", method = "lm")),
        title = "Pairwise Relationships in Iris Dataset")

Interpretation:

Off-diagonal panels: Reveal linear or non-linear associations. For instance, petal length and petal width show a strong positive correlation across species.

Diagonal panels: Show that Sepal.Width is roughly normally distributed, while Petal.Length exhibits bimodal/multimodal tendencies due to the separation of the setosa species from the others.

Group separation: Setosa (red) is clearly separable from the other two species across almost all variable combinations.

3. Correlation Heatmaps

A correlation matrix heatmap visualizes the correlation coefficients (\(r\)) between all pairs of numeric variables. Values range from -1 (perfect negative correlation) to +1 (perfect positive correlation).

# Compute correlation matrix for numeric variables
cor_matrix <- cor(iris[, 1:4])

# Melt the correlation matrix for ggplot2
melted_cor <- melt(cor_matrix)

# Plot heatmap
ggplot(melted_cor, aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-1, 1), name = "Pearson\nCorrelation") +
  theme_minimal() +
  labs(title = "Correlation Heatmap of Iris Features",
       x = "", y = "") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

Interpretation:

Red cells indicate strong positive correlations (e.g., Petal Length and Petal Width, \(r \approx 0.96\)).

Blue/White cells indicate weaker or negative relationships (e.g., Sepal Width has a weak negative correlation with Petal Length and Petal Width).

4. Parallel Coordinates Plots

Parallel coordinates map each multivariate observation to a polygonal line. Each variable corresponds to a vertical axis, and the values are scaled (usually between 0 and 1) so they fit on a uniform axis. This method excels at spotting clusters, outliers, and trends across many dimensions simultaneously.

# Load MASS library for parcoord
library(MASS)

# Parallel coordinates plot colored by species
parcoord(iris[, 1:4], col = as.numeric(iris$Species), 
         var.label = TRUE, 
         main = "Parallel Coordinates Plot of Iris Data")
legend("topright", legend = levels(iris$Species), col = 1:3, lty = 1, cex = 0.8)

Interpretation:

Trend profiles: Each colored line represents a single flower.

Feature discriminators: Setosa (black lines) clearly drops low on Petal Length and Petal Width axes compared to versicolor and virginica, indicating these dimensions are key differentiators.

5. Radar Charts (Spider Plots or Star Plots)

A radar chart (also known as a spider or star plot) projects multiple numerical metrics along independent axes extending outward from a central origin. Each observation forms a closed polygonal shape, making it intuitive to compare multi-attribute profiles against a baseline or across different categorical groups.

Because standard base R packages do not include a native radar chart function for data.frame objects, we typically use the fmsb package, which requires data to be formatted with specific minimum and maximum rows at the top.

# Install and load the fmsb package if not already installed
# install.packages("fmsb")
library(fmsb)

# 1. Prepare data: Calculate mean feature values for each Iris species
agg_iris <- aggregate(iris[, 1:4], by = list(Species = iris$Species), FUN = mean)

# 2. Extract numeric features and set row names to species
radar_data <- agg_iris[, 2:5]
rownames(radar_data) <- agg_iris$Species

# 3. fmsb requires maximum and minimum rows for the radar axes scaling
# Let's define bounds (min = 0, max = max value across all features + buffer)
max_vals <- apply(radar_data, 2, max)
min_vals <- apply(radar_data, 2, min)

# Combine bounds with the actual data
radar_df <- rbind(max_vals, min_vals, radar_data)

# 4. Define colors and aesthetics for plotting
colors_border <- c("#E41A1C", "#377EB8", "#4DAF4A")
colors_in <- c(scales::alpha("#E41A1C", 0.2),
               scales::alpha("#377EB8", 0.2),
               scales::alpha("#4DAF4A", 0.2))

# 5. Generate the radar chart
par(mar = c(1, 1, 2, 1))
radarchart(radar_df, 
           axistype = 1,
           # Customize polygon colors
           pcol = colors_border, 
           pfcol = colors_in, 
           plwd = 2, 
           plty = 1,
           # Customize grid aesthetics
           cglcol = "grey", 
           cglty = 1, 
           axislabcol = "grey", 
           cglwd = 0.8,
           # Labels
           title = "Radar Chart of Iris Species Mean Profiles")

# Add a legend
legend("topright", 
       legend = rownames(radar_data), 
       bty = "n", pch = 20, col = colors_border, 
       text.col = "black", pt.cex = 1.5)

Interpretation:

Polygon shape & area: Each species forms a distinct geometric footprint. For example, setosa (red polygon) forms a tiny, skewed shape restricted mostly to the sepal width axis, showing low values for petal lengths and widths.

Profile comparison: Virginica (green polygon) spans the largest overall surface area, pushing outward significantly along the Petal.Length and Petal.Width axes.

Multi-attribute evaluation: Unlike scatterplots that map relationships between pairs, the radar chart displays all four variables simultaneously for group profiling, making it ideal for stakeholder presentations where quick visual comparisons of overall performance or traits are needed.

6. Bubble Charts

A bubble chart extends a standard scatter plot by utilizing position (X and Y axes), marker size (representing a 3rd quantitative variable), and color shading (representing a 4th categorical or continuous variable).

ggplot(iris, aes(x = Petal.Length, y = Petal.Width, size = Sepal.Width, color = Species)) +
  geom_point(alpha = 0.6) +
  scale_size_continuous(range = c(1, 6)) +
  theme_minimal() +
  labs(title = "Multivariate Bubble Chart",
       subtitle = "X: Petal.L | Y: Petal.W | Size: Sepal.W | Color: Species",
       x = "Petal Length", y = "Petal Width", size = "Sepal Width")

Interpretation:

Bubble chart allows us to visualize interactions across 4 dimensions simultaneously, illustrating that larger petals (versicolor and virginica) also tend to fluctuate more noticeably in sepal width size gradients.

7. Principal Component Analysis (PCA) Biplot

PCA is a dimensionality reduction technique that transforms correlated multivariate data into a smaller set of uncorrelated variables called principal components (PCs). A biplot displays both the observations (as points) and the original variables (as vectors) projected onto the first two principal components.

# Perform PCA with standard scaling
pca_res <- prcomp(iris[, 1:4], scale. = TRUE)

# Calculate percentage of variance explained per component
var_explained <- round(pca_res$sdev^2 / sum(pca_res$sdev^2) * 100, 1)

# Compile projection frame
pca_df <- data.frame(pca_res$x, Species = iris$Species)

# Render PCA scatter with confidence ellipses
ggplot(pca_df, aes(x = PC1, y = PC2, color = Species)) +
  geom_point(size = 2, alpha = 0.8) +
  stat_ellipse(level = 0.95) +
  labs(title = "PCA Biplot: Low-Dimensional Projection",
       x = paste0("PC1 (", var_explained[1], "%)"),
       y = paste0("PC2 (", var_explained[2], "%)")) +
  theme_minimal()

Interpretation:

Variance captured: PC1 and PC2 typically capture over 95% of the total variance in the iris dataset, allowing us to visualize 4D data effectively in 2D.

Cluster separation: Setosa forms a distinct cluster on the left, while versicolor and virginica overlap slightly along PC2.

Ellipses: The 95% confidence concentration ellipses confirm tight grouping within species and clear separation between setosa and the rest.