Hotelling’s T-Squared Ellipse
Calculating the semi-axes
To visualize the confidence region for our multivariate data, we
employ the ellipseParam() function to generate a confidence
ellipse. Our objective is to calculate the lengths of the semi-axes for
this ellipse, focusing on the bivariate relationship within the PC1-PC3
subspace of our principal component analysis. We maintain the default
value for k (number of components) at 2. This ensures we’re
working with a two-dimensional representation. We specify
pcx = 1 and pcy = 3 as inputs. This directs
the function to use the 1st and 3rd principal components for the
x and y axes, respectively. The Hotelling’s T-squared
statistic returned in Tsquare is computed on these same two
components, so an observation lies outside the ellipse exactly when its
T-squared value exceeds the corresponding cutoff.
str(ellipse_axes)
#> List of 5
#> $ Tsquare : tibble [100 × 1] (S3: tbl_df/tbl/data.frame)
#> ..$ value: num [1:100] 17.124 1.195 0.818 2.286 0.392 ...
#> $ cutoff.99pct: num 9.76
#> $ cutoff.95pct: num 6.24
#> $ nb.comp : int 2
#> $ Ellipse : tibble [1 × 5] (S3: tbl_df/tbl/data.frame)
#> ..$ a.99pct: num 19369
#> ..$ b.99pct: num 8416
#> ..$ a.95pct: num 15492
#> ..$ b.95pct: num 6732
#> ..$ angle : num 0We can extract parameters for further use:
- Semi-axes of the ellipse at 99% confidence level.
a1 <- ellipse_axes %>% pluck("Ellipse", "a.99pct")
b1 <- ellipse_axes %>% pluck("Ellipse", "b.99pct")- Semi-axes of the ellipse at 95% confidence level.
a2 <- ellipse_axes %>% pluck("Ellipse", "a.95pct")
b2 <- ellipse_axes %>% pluck("Ellipse", "b.95pct")- Hotelling’s T\(^2\) statistic.
t1 <- round(as.numeric(pca_mod$eig[1,2]), 2)
t2 <- round(as.numeric(pca_mod$eig[2,2]), 2)
t3 <- round(as.numeric(pca_mod$eig[3,2]), 2)pca_scores %>%
ggplot(aes(x = Dim.1, y = Dim.3)) +
geom_ellipse(aes(x0 = 0, y0 = 0, a = a1, b = b1, angle = 0), linewidth = .5, linetype = "solid", fill = "white") +
geom_ellipse(aes(x0 = 0, y0 = 0, a = a2, b = b2, angle = 0), linewidth = .5, linetype = "solid", fill = "white") +
geom_point(aes(fill = Tsq), shape = 21, size = 3, color = "black") +
scale_fill_viridis_c(option = "viridis") +
geom_hline(yintercept = 0, linetype = "solid", color = "black", linewidth = .2) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = .2) +
labs(
title = "Scatterplot of PCA scores",
subtitle = "PC1 vs. PC3",
x = glue("PC1 [{t1}%]"),
y = glue("PC3 [{t3}%]"),
fill = "T2"
) +
theme_grey() +
theme(
aspect.ratio = .7,
panel.grid = element_blank(),
panel.background = element_rect(
colour = "black",
linewidth = .3
)
)
Calculating the x-y coordinates
An alternative method for visualizing the Hotelling’s T-squared
confidence region is to utilize the ellipseCoord()
function. This function generates the x and y
coordinates for plotting the confidence ellipse, offering greater
flexibility in visualization and further analysis. It allows users to
specify a custom confidence level via the confi.limit
parameter. The default confidence level is set at 95%, which is commonly
used in statistical analyses. The function returns a set of coordinates
that define the ellipse’s boundary in the chosen subspace, thereby
complementing the semi-axes information provided by the
ellipseParam() function.
In the example below, we focus on the subspace spanned by the 2nd and 3rd components.
str(xy_coord)
#> tibble [500 × 2] (S3: tbl_df/tbl/data.frame)
#> $ x: num [1:500] 9620 9620 9617 9613 9608 ...
#> $ y: num [1:500] -3.22e-12 9.44e+01 1.89e+02 2.83e+02 3.77e+02 ...ggplot() +
geom_polygon(data = xy_coord, aes(x, y), color = "black", fill = "white") +
geom_point(data = pca_scores, aes(x = Dim.2, y = Dim.3), shape = 21, size = 3, fill = "black", color = "black", alpha = 0.7) +
geom_hline(yintercept = 0, linetype = "solid", color = "black", linewidth = .2) +
geom_vline(xintercept = 0, linetype = "solid", color = "black", linewidth = .2) +
labs(
title = "Scatterplot of PCA scores",
subtitle = "PC2 vs. PC3",
x = glue("PC2 [{t2}%]"),
y = glue("PC3 [{t3}%]")
) +
theme_grey() +
theme(
aspect.ratio = .7,
panel.grid = element_blank(),
panel.background = element_rect(
colour = "black",
linewidth = .3
)
)