data(USArrests)
X <- USArrests                        # Murder, Assault, UrbanPop, Rape (all numeric)

res.pca <- prcomp(X, center = TRUE, scale = TRUE)   # PCA via SVD, z-standardised internally

# Sign convention (RAISINS default): reverse sign for consistent orientation
res.pca$x        <- -res.pca$x
res.pca$rotation <- -res.pca$rotation

res.pca$x                              # PC scores (individuals)
res.pca$rotation                       # loadings (variables)
res.pca$sdev                           # standard deviations of the PCs

library(factoextra)
library(scales)

# Eigenvalues and % variance explained
eig.val <- get_eigenvalue(res.pca)
rownames(eig.val) <- paste0("PC", seq_len(nrow(eig.val)))
colnames(eig.val) <- c("Eigen value", "Variance(%)", "Cumulative variance")

e1 <- eig.val[1, "Variance(%)"]        # % variance explained by PC1
e2 <- eig.val[2, "Variance(%)"]        # % variance explained by PC2

# Loadings (variable contributions to each PC)
load <- as.data.frame(res.pca$rotation)

# Index score - PC1 and PC2 scores rescaled to [0, 1]
index1 <- res.pca$x[, 1]
index2 <- res.pca$x[, 2]

index1_scaled <- rescale(index1, to = c(0, 1))
index2_scaled <- rescale(index2, to = c(0, 1))

index_1 <- data.frame(`Index Score` = index1, `Scaled Index` = index1_scaled, check.names = FALSE)
index_2 <- data.frame(`Index Score` = index2, `Scaled Index` = index2_scaled, check.names = FALSE)

library(cluster)

# Use the first 3 PCs (as displayed on the 2D biplot / 3D scatter)
df <- as.data.frame(res.pca$x[, 1:3])

compute_sil <- function(mat) {
  kmax <- min(10, nrow(mat) - 1)
  sil  <- rep(NA, kmax)
  for (k in 2:kmax) {
    km <- kmeans(mat, centers = k)
    s  <- cluster::silhouette(km$cluster, dist(mat))
    sil[k] <- mean(s[, 3])
  }
  which.max(sil)
}

k12 <- compute_sil(df[, c(1, 2)])      # PC1 vs PC2 (2D biplot default)
k13 <- compute_sil(df[, c(1, 3)])      # PC1 vs PC3
k23 <- compute_sil(df[, c(2, 3)])      # PC2 vs PC3

suggested_k <- na.omit(c(k12, k13, k23))[1]

# Fallback: WSS elbow on (PC1, PC2) if silhouette evaluation fails
if (length(suggested_k) == 0) {
  wss <- sapply(1:5, function(k) sum(kmeans(df[, 1:2], centers = k)$withinss))
  suggested_k <- which.min(diff(wss))
}

# Final fallback
if (length(suggested_k) == 0 || is.na(suggested_k) || suggested_k < 2) suggested_k <- 2

# Cluster assignment used on the 2D biplot
km_biplot <- kmeans(df[, c("PC1", "PC2")], centers = suggested_k)

library(factoextra)

pca_var <- get_pca_var(res.pca)        # variable coordinates, contributions, cos2
pca_ind <- get_pca_ind(res.pca)        # individual coordinates, cos2

df_var <- as.data.frame(pca_var$coord[, 1:2])          # variable arrow endpoints (PC1, PC2)
df_var$contrib <- rowSums(pca_var$contrib[, 1:2])       # combined contribution to PC1+PC2
df_var$cos2    <- rowSums(pca_var$cos2[, 1:2])           # combined quality of representation

df_ind <- as.data.frame(pca_ind$coord[, 1:2])            # individual (treatment) coordinates
df_ind$cos2 <- rowSums(pca_ind$cos2[, 1:2])
df_ind$cluster <- factor(km_biplot$cluster)               # cluster overlay from previous step

# Correlation circle uses the same variable coordinates (PC1, PC2),
# plotted against a unit reference circle
circle_var <- as.data.frame(pca_var$coord[, 1:2])