Skip to content

Commit

Permalink
Minor updates
Browse files Browse the repository at this point in the history
  • Loading branch information
ChristianGoueguel committed Apr 25, 2024
1 parent a9ab0fb commit ae30e2a
Show file tree
Hide file tree
Showing 2 changed files with 5 additions and 38 deletions.
14 changes: 0 additions & 14 deletions R/ellipseCoord.R
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,6 @@
#' xy_coord <- ellipseCoord(data = pca_scores, pcx = 1, pcy = 2, conf.limit = 0.95, pts = 200)
#'
ellipseCoord <- function(data, pcx = 1, pcy = 2, conf.limit = 0.95, pts = 200) {

# check input validity
if (length(data) == 0) {
stop("Data must not be empty.")
}
Expand All @@ -44,32 +42,20 @@ ellipseCoord <- function(data, pcx = 1, pcy = 2, conf.limit = 0.95, pts = 200) {
if (conf.limit < 0 || conf.limit > 1) {
stop("Confidence level should be between 0 and 1.")
}

# matrix of data
X <- as.matrix(data)

# Sample size
n <- nrow(X)

# Confidence limit
alpha <- as.numeric(conf.limit)

# Number of points
m <- as.numeric(pts)
p <- seq(0, 2*pi, length.out = m)

# # Hotelling’s T-square limit
Tsq_limit <- (2*(n-1)/(n-2))*stats::qf(p = alpha, df1 = 2, df2 = (n-2))

# Coordinate points
rx <- sqrt(Tsq_limit*stats::var(X[, pcx]))
ry <- sqrt(Tsq_limit*stats::var(X[, pcy]))

res.coord <- tibble::tibble(
x = rx*cos(p) + mean(X[, pcx], na.rm = TRUE),
y = ry*sin(p) + mean(X[, pcy], na.rm = TRUE)
)

return(res.coord)

}
29 changes: 5 additions & 24 deletions R/ellipseParam.R
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
#'4. **cutoff.95pct** Number corresponding to the T-squared cutoff at 95% confidence level.
#' @export ellipseParam
#' @author Christian L. Goueguel,
#' christian.goueguel@gmail.com
#' @examples
#' ## Principal components analysis (PCA)
#' library(dplyr)
Expand All @@ -31,8 +30,6 @@
#' T2 <- ellipseParam(data = pca_scores, k = 2, pcx = 1, pcy = 2)
#'
ellipseParam <- function(data, k = 2, pcx = 1, pcy = 2) {

# check input validity
if (length(data) == 0) {
stop("Data must not be empty.")
}
Expand All @@ -45,52 +42,37 @@ ellipseParam <- function(data, k = 2, pcx = 1, pcy = 2) {
if (pcx == 0 || pcy == 0) {
stop("pcx and pcy must be non-zero.")
}
if (k < 2) {
stop("k must be at least equal to 2.")
if (!is.numeric(k) || k < 2 || k != floor(k)) {
stop("k must be a positive integer >= 2.")
}
if (k > ncol(data)) {
stop("k exceeds the number of component in the data.")
stop("k exceeds the number of columns in the data.")
}

# matrix of data
X <- as.matrix(data)

# sample size
n <- nrow(X)

# number of principal component
A <- as.numeric(k)

# Squared Mahalanobis distance
MDsq <- stats::mahalanobis(
x = X,
center = colMeans(X),
cov = stats::cov(X),
inverted = FALSE
)

if(k > 2) {
# 99% and 95% confidence limit
Tsq_limit1 <- (A*(n-1)/(n-A))*stats::qf(p = 0.99, df1 = A, df2 = (n-A))
Tsq_limit2 <- (A*(n-1)/(n-A))*stats::qf(p = 0.95, df1 = A, df2 = (n-A))

# Hotelling’s T-square
Tsq <- tibble::tibble(value = ((n-A)/(A*(n-1)))*MDsq)

res_list <- list(
"Tsquare" = Tsq,
"cutoff.99pct" = Tsq_limit1,
"cutoff.95pct" = Tsq_limit2
)

return(res_list)
}

if(k == 2) {
# 99% and 95% confidence limit
Tsq_limit1 <- (A*(n-1)/(n-A))*stats::qf(p = 0.99, df1 = A, df2 = (n-A))
Tsq_limit2 <- (A*(n-1)/(n-A))*stats::qf(p = 0.95, df1 = A, df2 = (n-A))

# Hotelling ellipse semi-axes
a_limit1 <- sqrt(Tsq_limit1*stats::var(X[, pcx]))
a_limit2 <- sqrt(Tsq_limit2*stats::var(X[, pcx]))
b_limit1 <- sqrt(Tsq_limit1*stats::var(X[, pcy]))
Expand All @@ -103,7 +85,6 @@ ellipseParam <- function(data, k = 2, pcx = 1, pcy = 2) {
b.95pct = b_limit2
)

# Hotelling’s T-square
Tsq <- tibble::tibble(value = ((n-2)/(2*(n-1)))*MDsq)

res_list <- list(
Expand All @@ -112,7 +93,7 @@ ellipseParam <- function(data, k = 2, pcx = 1, pcy = 2) {
"cutoff.99pct" = Tsq_limit1,
"cutoff.95pct" = Tsq_limit2
)

return(res_list)
}
}

0 comments on commit ae30e2a

Please sign in to comment.