## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  out.width = "75%",
  fig.align = "center",
  fig.width = 8,
  fig.height = 5,
  dpi = 96,
  comment = "#>"
)

## ----setup, include=FALSE-----------------------------------------------------
library(ggplot2)
library(gghotelling)
library(MASS)

theme_set(theme_minimal())

## ----example------------------------------------------------------------------
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()

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

## ----example2-----------------------------------------------------------------
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"))

## ----example2b----------------------------------------------------------------
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)))

## ----example3-----------------------------------------------------------------
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)


## ----example3b----------------------------------------------------------------
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)

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

head(outlier_stats)

## ----outlier_plot,fig.width=8, fig.height=4-----------------------------------
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)")

## ----outlier_plot-2,eval=FALSE------------------------------------------------
# plot_outliers(pca_df[ , c("PC1", "PC2")], level = 0.95)

## -----------------------------------------------------------------------------
library(HDclassif)
data(wine)
wine <- wine[ wine$class == 1, ]
wine <- data.frame("malic_acid"=wine$V2, "proline"=wine$V13)

ggplot(wine, aes(malic_acid, proline)) +
  geom_hotelling(type="c2data", level = .975, color = "red") +
  geom_hotelling(type="c2data", level = .975, robust = TRUE, color = "blue") +
  geom_point() +
  annotate("text", x=2.5, y = 1675, label = "MCD", color = "blue") +
  annotate("text", x=3.5, y = 1400, label = "Classical", color = "red")

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

## -----------------------------------------------------------------------------
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()

## -----------------------------------------------------------------------------
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)

## -----------------------------------------------------------------------------
# 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()

## ----example_autoplot---------------------------------------------------------
autoplot(pca, group = iris$Species) + 
  autolayer(pca, group = iris$Species)

## -----------------------------------------------------------------------------
pca <- prcomp(iris[, 1:4], scale.=TRUE)
plot(pca$x[, 1:2], col=c(3,4,6)[iris$Species], pch=19)

# create an ellipse for each species
eli <- tapply(1:nrow(iris), iris$Species,
              \(vec) hotelling_ellipse(pca$x[vec,1:2], level=.95))

# plot the ellipse
lapply(1:3, \(i) lines(eli[[i]], col=c(3,4,6)[i], lwd=2))

# create outlier stats for each point in each species
outl <- tapply(1:nrow(iris), iris$Species,
              \(vec) outliers(pca$x[vec,1:2], level=.95))

# make points and labels
outl <- cbind(pca$x[, 1:2], do.call(rbind, outl))
points(outl[outl$is_outlier, 1:2], col="red", pch=19, cex=2)
text(outl[outl$is_outlier, 1:2], labels=rownames(outl)[outl$is_outlier],
     pos=3, cex=.8)