Installation

You can install the development version of gghotelling from GitHub with:

# install.packages("pak")
pak::pak("january3/gghotelling")

# or devtools::install_github("january3/gghotelling")

The geom_hotelling() function

The package defines a new geom, geom_hotelling(), which can be used to add Hotelling ellipses to ggplot2 scatter plots.

library(ggplot2)
library(gghotelling)

pca <- prcomp(iris[, 1:4], scale.=TRUE)
pca_df <- cbind(iris, pca$x)

# set the confidence level
ggplot(pca_df, aes(PC1, PC2)) +
  geom_hotelling(level=.99) +
  geom_point()

ggplot(pca_df, aes(PC1, PC2, color=Species)) +
  geom_hotelling() +
  geom_point()

ggplot(pca_df, aes(PC1, PC2, color=Species)) +
  geom_hotelling(alpha=0.1, aes(fill = Species)) +
  geom_point()

# set custom CI/coverage level
ggplot(pca_df, aes(PC1, PC2, color=Species)) +
  geom_hotelling(alpha=0.1, aes(fill = Species), level=.99) +
  geom_point()

Types of Ellipses

How are the Hotelling ellipses different, say, from the ellipses created by stat_ellipse() in ggplot2, or the ellipse::ellipse() function?

Actually, the geom_hotelling() function can create three different types of ellipses with two flavors:

• Hotelling T² data ellipses (default, type="t2data"): these ellipses represent the spread of the data points themselves, based on the Hotelling T² distribution. They can be used to visualize the overall distribution of the data and identify potential outliers.
• Hotelling T² confidence ellipses for group means (type="t2mean"): these ellipses represent the confidence region for the mean of each group, based on the Hotelling T² distribution. They can be used to compare the means of different groups and assess whether they are significantly different from each other.
• Chi-squared data ellipses (type="c2data"): these ellipses are based on the chi-squared distribution and also represent the spread of the data points. They are similar to the ellipses created by stat_ellipse() in ggplot2 and the ellipse::ellipse() function.
• Each of the above can be produced as a robust version using the MCD estimator by setting robust=TRUE.

All three ellipses above use Mahalanobis distance contours, but differ in the statistical choice of distribution (Hotelling T² vs χ²) in order to select the Mahalanobis distance threshold for drawing the ellipse.

So why the different distributions? The point is whether we are considering the data to be the actual population (in which case we use the χ² distribution) or a sample from a larger population (in which case we use the Hotelling T² distribution). The Hotelling T² distribution takes into account the uncertainty in estimating the population parameters (mean and covariance) from a finite sample, leading to wider ellipses compared to the χ² distribution, as you can see on the figure below - the dashed ellipses are the χ² data ellipses:

ggplot(pca_df, aes(PC1, PC2, color = Species)) +
  geom_hotelling(level=.99) +
  geom_hotelling(level=.99, type="c2data", linetype = "dashed") +
  geom_point()

In addition to the classical Hotelling ellipses, robust versions can be created with the robust=TRUE argument, which uses the Minimum Covariance Determinant (MCD) estimator to compute robust estimates of the mean and covariance matrix (see below for details).

Outlier Detection

The package also provides per-point, group-wise T² statistics which can be used to identify multivariate outliers.

ggplot(pca_df, aes(PC1, PC2, group=Species)) +
  geom_hotelling(level = 0.75, alpha=0.1, aes(fill = Species)) +
  # add points and calculate outlier stats; we assign the `is_outlier` variable
  # calculated by stat_outliers() to the color aesthetic
  stat_outliers(level = .75, 
                aes(shape = Species, color = after_stat(is_outlier))) +
  # color outliers in red
  scale_color_manual(values=c("TRUE"="red", "FALSE"="grey"))

The stat_outliers() function computes the Hotelling T² statistic for each point and compares it to the critical value from the Hotelling T² distribution at the specified confidence level (default is 0.95). The new data is available to aesthetics via after_stat(), including:

• d2: the squared Mahalanobis distance for each point
• c2crit: the critical value for d2 from the Chi-squared distribution at the specified confidence level
• t2crit: the critical value for d2 from the Hotelling T2 distribution at the specified confidence level
• is_outlier: a logical indicating whether the point is an outlier

These variables can be used, through stat_outliers(), to map aesthetics such as color, shape, or size to highlight outliers. For example:

ggplot(pca_df, aes(PC1, PC2, group=Species)) +
  geom_hotelling(alpha=0.1, aes(fill = Species)) +
  stat_outliers(size=2, aes(shape = Species, color = after_stat(d2)))

This can be useful for identifying potential outliers in multivariate data. The outliers can be directly labeled as follows:

ggplot(pca_df, aes(PC1, PC2, group=Species, label=rownames(pca_df))) +
  geom_hotelling(alpha=0.1, aes(fill = Species)) +
  geom_point(aes(color = Species)) +
  stat_outliers(geom="label", 
                        outlier_only = TRUE)

Or even better, using ggrepel to avoid overlapping labels:

library(ggrepel)
ggplot(pca_df, aes(PC1, PC2, group=Species, label=rownames(pca_df))) +
  geom_hotelling(alpha=0.1, aes(fill = Species)) +
  geom_point(aes(color = Species)) +
  stat_outliers(geom="label_repel",
                        outlier_only = TRUE)

The actual calculation of the Hotelling T² statistics and critical values is done in the function outliers(), which can also be used directly on data frames to compute the statistics without plotting:

outlier_stats <- outliers(pca_df[ , c("PC1", "PC2")], level = 0.95)

head(outlier_stats)
#>         d2   t2crit   c2crit is_outlier
#> 1 1.996071 6.155707 5.991465      FALSE
#> 2 1.967770 6.155707 5.991465      FALSE
#> 3 2.029500 6.155707 5.991465      FALSE
#> 4 2.187372 6.155707 5.991465      FALSE
#> 5 2.398597 6.155707 5.991465      FALSE
#> 6 3.876401 6.155707 5.991465      FALSE

We can visualize it with the typical ggplot2 syntax:

outlier_stats$id <- 1:nrow(outlier_stats)
outlier_labels <- ifelse(outlier_stats$is_outlier,
                             as.character(outlier_stats$id), NA)
ggplot(outlier_stats, aes(x = id, y = sqrt(d2))) +
  geom_segment(aes(xend = id, yend = 0), alpha = .3) +
  geom_point(aes(color = is_outlier), size = 2) +
  scale_color_manual(values=c("TRUE"="red", "FALSE"="black")) +
  geom_label(aes(label = outlier_labels), nudge_y = 0.2, na.rm = TRUE) +
  geom_hline(aes(yintercept = sqrt(t2crit)), color = "red", linetype = "dashed") +
  annotate("text", x = 1, y = sqrt(outlier_stats$t2crit[1]) + 0.1,
           label = "Critical value", color = "red", hjust = 0) +
  theme(legend.position = "none") +
  labs(y = "Mahalanobis distance (T² statistic)")

For convenience, there is a plot_outliers() function that creates the above plot directly from a data frame:

plot_outliers(pca_df[ , c("PC1", "PC2")], level = 0.95)

Convex hulls with geom_hull()

The package provides basic convex hull:

ggplot(iris, aes(Sepal.Length, Sepal.Width, color=Species)) +
  geom_hull(mapping = aes(fill = Species), alpha=.1) +
  geom_point()

Bagplots with geom_bag()

With geom_bag() we can create bagplots, which are bivariate generalizations of boxplots. The bagplot shows the central “bag” containing 50% of the data points, the “loop” which is an expanded region that helps identify outliers, and the outliers themselves. The geom_bag() can plot either the bag or the loop; plotting both together requires calling geom_bag() twice:

ggplot(iris, aes(Sepal.Width, Sepal.Length, color=Species)) +
  geom_bag(aes(fill=Species), alpha=.3) +
  geom_bag(aes(fill=Species), alpha=.1, what = "loop") +
  geom_point()

Kernel density estimate contours with geom_kde()

Furthermore, we also have geom_kde() to draw 2D kernel density estimate. Unlike geom_density_2d(), this function actually creates a single contour for the specified coverage; for example, if coverage=0.95, the contour encloses roughly 95% of the data points. This is useful for visualizing the density distribution of non-elliptical data:

df <- data.frame(x=rnorm(500) + 5)
df$y <- df$x^5 + rnorm(500)*500
ggplot(df, aes(x=x, y=y)) +
  geom_point()+
  geom_hull(color = "grey") +
  geom_kde(color="red", linewidth=1)

As you can see, the geom_kde() contour nicely captures the non-elliptical distribution of the data, while the convex hull includes a large empty area.

Of course, you can also use aesthetics like fill with geom_kde() and overlay several contours:

# interesting little fact: ggplot2 happily accepts lists of geoms/layers
# and adds them one by one to the plot

rings <- lapply(seq(.05, .95, length.out = 10), \(i) {
  geom_kde(aes(fill = Species), alpha =.05, coverage = i)
})

ggplot(pca_df, aes(x = PC1, y = PC2, color=Species)) +
  rings +
  geom_point()