This function calculates Hotelling’s T-squared statistic and, when applicable, the lengths of the semi-axes of the Hotelling’s ellipse. It can work with a specified number of components or use a cumulative variance threshold.
A matrix, data frame or tibble containing scores from PCA, PLS, ICA, or other similar methods. Each column should represent a component, and each row an observation.
An integer specifying the number of components to use (default is 2). This parameter is ignored if threshold is provided.
An integer specifying which component to use for the x-axis when k = 2 (default is 1).
An integer specifying which component to use for the y-axis when k = 2 (default is 2).
A numeric value between 0 and 1 specifying the desired cumulative explained variance threshold (default is NULL). If provided, the function determines the minimum number of components needed to explain at least this proportion of total variance. When NULL, the function uses the fixed number of components specified by k.
A numeric value specifying the minimum proportion of total variance a component should explain to be considered non-negligible (default is 0.001, i.e., 0.1%).
A numeric value specifying the minimum absolute variance a component should have to be considered non-negligible (default is .Machine$double.eps).
A character string specifying how the T-squared cutoffs are computed: "f" (default) or "beta". See Details.
A numeric vector of confidence levels, each strictly between 0 and 1, at which the T-squared cutoffs and ellipse semi-axes are computed (default is c(0.95, 0.99)). Any number of levels can be given. See Value for how the results are named.
A list containing the following elements:
Tsquare: A data frame containing the T-squared statistic for each observation (the squared Mahalanobis distance), on the same scale as the cutoffs. When k = 2, it is computed on components pcx and pcy.
cutoff.<level>pct: The T-squared cutoff value at each confidence level in conf.limit, from the highest to the lowest level. With the default conf.limit = c(0.95, 0.99), these are cutoff.99pct and cutoff.95pct; with conf.limit = c(0.975, 0.999), they are cutoff.99.9pct and cutoff.97.5pct.
nb.comp: The number of components used in the calculation.
Ellipse: A data frame containing the lengths of the semi-axes at each confidence level (a.<level>pct, b.<level>pct, e.g. a.99pct, b.99pct, a.95pct, b.95pct by default) and the rotation angle of the ellipse in radians (only when k = 2). a is the semi-axis closest to the pcx direction. For uncorrelated scores, such as the PCA or PLS scores of the samples the model was fitted on, angle is 0.
When threshold is used, the function selects the minimum number of k components
that cumulatively explain at least the specified proportion of variance. This
parameter allows for dynamic component selection based on explained variance,
rather than using a fixed number of components. It must be greater than rel.tol.
Typical values range from 0.8 to 0.95.
The rel.tol parameter sets a minimum variance threshold for individual components.
Components with variance below this threshold are considered negligible and are
removed from the analysis. Setting rel.tol too high
may remove potentially important components, while setting it too low may
retain noise or cause computational issues. Adjust based on your data
characteristics and analysis goals.
Note that components are considered to have near-zero variance and are removed
if their relative variance is below rel_tol or their absolute variance is
below abs_tol. This dual-threshold approach helps ensure numerical stability
while also accounting for the relative importance of components. The default
value for abs.tol is set to .Machine$double.eps, providing a lower bound
for detecting near-zero variance that may cause numerical instability.
The method parameter sets the distribution used for the cutoffs. With
method = "f" (default), the cutoffs are
\(\frac{k(n-1)}{n-k} F_{\alpha}(k, n-k)\),
as in previous versions of the package. With method = "beta", the cutoffs
follow the exact distribution of T-squared for the observations used to
estimate the mean and covariance,
\(\frac{(n-1)^2}{n} B_{\alpha}(k/2, (n-k-1)/2)\)
(Tracy, Young and Mason, 1992), e.g. the scores of the samples a PCA or PLS
model was built on. The F-based limit is more conservative, especially for
small n: it can even exceed the largest T-squared value any observation can
reach, \((n-1)^2/n\). For n larger than about 100, the two
limits are close.
When the selected components are correlated (e.g. new samples projected onto a model, or ICA scores), the ellipse is rotated so that it matches the T-squared statistic: an observation lies outside the ellipse exactly when its T-squared value exceeds the cutoff.
Tracy, N. D., Young, J. C. and Mason, R. L. (1992). Multivariate control charts for individual observations. Journal of Quality Technology, 24(2), 88–95.
if (FALSE) { # \dontrun{
# Load required libraries
library(HotellingEllipse)
library(dplyr)
data("specData", package = "HotellingEllipse")
# Perform PCA
set.seed(123)
pca_mod <- specData %>%
select(where(is.numeric)) %>%
FactoMineR::PCA(scale.unit = FALSE, graph = FALSE)
# Extract PCA scores
pca_scores <- pca_mod$ind$coord %>% as.data.frame()
# Example 1: Calculate Hotelling’s T-squared and ellipse parameters using
# the 2nd and 4th components
T2_fixed <- ellipseParam(x = pca_scores, pcx = 2, pcy = 4)
# Example 2: Calculate using the first 4 components
T2_comp <- ellipseParam(x = pca_scores, k = 4)
# Example 3: Calculate using a cumulative variance threshold
T2_threshold <- ellipseParam(x = pca_scores, threshold = 0.95)
# Example 4: Cutoffs and semi-axes at the 97.5% and 99.9% confidence levels
T2_levels <- ellipseParam(x = pca_scores, conf.limit = c(0.975, 0.999))
T2_levels$cutoff.99.9pct
T2_levels$Ellipse$a.97.5pct
} # }