Implementing Horn's Parallel Analysis for Covariance PCA

# Load necessary libraries
library(nFactors)
library(psych)

# Import your data
data <- read.xlsx("data.xlsx")

# Ensure your data is properly scaled
# Important: Do NOT standardize the data for covariance matrix PCA
n_subjects <- nrow(data)
n_variables <- ncol(data)

# Calculate the covariance matrix
cov_mat <- cov(data, use = "pairwise.complete.obs")

# Calculate eigenvalues from your actual data
eigen_observed <- eigen(cov_mat)$values

# Perform parallel analysis for covariance matrix
results <- parallel(
  var = n_variables,
  subject = n_subjects,
  rep = 1000,          # Number of replications
  quantile = 0.95,     # 95th percentile cutoff
  model = "components", # PCA model
  sd = sqrt(diag(cov_mat)) # Critical: use actual standard deviations
)

# Extract simulated eigenvalues
simulated_eigen_mean <- results$eigen$mevpea
simulated_eigen_quantile <- results$eigen$qevpea

# Create comparison dataframe
comparison <- data.frame(
  Component = 1:n_variables,
  Observed = eigen_observed,
  Simulated_Mean = simulated_eigen_mean,
  Simulated_95th = simulated_eigen_quantile,
  Retain = eigen_observed > simulated_eigen_quantile
)

# View results
print(comparison)

library(paran)

# Perform parallel analysis
pa_result <- paran(
  cov_mat,
  n.iter = 1000,
  cent = 0.05,
  quantile = TRUE
)

# Plot results
plot(pa_result)

# Get suggested number of components
pa_result$nfact

# Manual parallel implementation
set.seed(123)  # For reproducibility
n_sim <- 1000
sim_eigen <- matrix(NA, nrow = n_variables, ncol = n_sim)

for (i in 1:n_sim) {
  # Generate random data with same covariance structure
  random_data <- matrix(rnorm(n_subjects * n_variables), 
                       nrow = n_subjects, ncol = n_variables)
  
  # Scale to match original variances (critical for covariance matrix)
  random_data_scaled <- sweep(random_data, 2, 
                             sqrt(diag(cov_mat)), `*`)
  
  # Calculate covariance matrix and eigenvalues
  cov_sim <- cov(random_data_scaled)
  sim_eigen[, i] <- eigen(cov_sim)$values
}

# Calculate mean and quantiles
sim_mean <- apply(sim_eigen, 1, mean)
sim_95th <- apply(sim_eigen, 1, quantile, probs = 0.95)

# Compare with observed eigenvalues
comparison <- data.frame(
  Component = 1:n_variables,
  Observed = eigen_observed,
  Simulated_Mean = sim_mean,
  Simulated_95th = sim_95th,
  Retain = eigen_observed > sim_95th
)

# Complete Horn's Parallel Analysis for Covariance Matrix PCA

# Load libraries
library(nFactors)
library(psych)
library(openxlsx)

# Set seed for reproducibility
set.seed(123)

# Import data
data <- read.xlsx("data.xlsx")

# Basic information about the data
cat("Number of observations:", nrow(data), "\n")
cat("Number of variables:", ncol(data), "\n")

# Calculate covariance matrix (do NOT standardize for covariance PCA)
cov_mat <- cov(data, use = "pairwise.complete.obs")

# Calculate observed eigenvalues
eigen_observed <- eigen(cov_mat)$values

# Perform parallel analysis
results <- parallel(
  var = ncol(data),
  subject = nrow(data),
  rep = 2000,          # Increased replications for stability
  quantile = 0.95,     # 95th percentile
  model = "components",
  sd = sqrt(diag(cov_mat))  # Critical parameter for covariance matrices
)

# Extract simulated eigenvalues
sim_mean <- results$eigen$mevpea
sim_95th <- results$eigen$qevpea

# Create comparison dataframe
comparison <- data.frame(
  Component = 1:length(eigen_observed),
  Observed_Eigenvalue = eigen_observed,
  Simulated_Mean = sim_mean,
  Simulated_95th = sim_95th,
  Difference = eigen_observed - sim_95th,
  Retain_Component = eigen_observed > sim_95th
)

# Display results
print("Eigenvalue Comparison:")
print(comparison)

# Determine number of components to retain
n_components <- sum(eigen_observed > sim_95th)
cat("\nRecommended number of components to retain:", n_components, "\n")

# Create visualization
plot(1:length(eigen_observed), eigen_observed, 
     type = "b", pch = 19, col = "blue",
     xlab = "Component Number", 
     ylab = "Eigenvalue",
     main = "Horn's Parallel Analysis - Covariance Matrix")
lines(1:length(eigen_observed), sim_mean, 
      type = "b", pch = 17, col = "red")
lines(1:length(eigen_observed), sim_95th, 
      type = "b", pch = 15, col = "green")
legend("topright", 
       legend = c("Observed", "Simulated Mean", "Simulated 95th"),
       col = c("blue", "red", "green"),
       pch = c(19, 17, 15))

# Additional diagnostic plots
par(mfrow = c(2, 2))
plotuScree(eigen_observed, main = "Scree Plot")
plotuScree(sim_95th, main = "Simulated 95th Percentile")
hist(eigen_observed, main = "Distribution of Observed Eigenvalues")
hist(sim_95th, main = "Distribution of Simulated 95th Percentile")
par(mfrow = c(1, 1))

Component 1: Observed = 5.2, Simulated_95th = 1.8 → RETAIN
Component 2: Observed = 2.1, Simulated_95th = 1.5 → RETAIN
Component 3: Observed = 1.2, Simulated_95th = 1.3 → DO NOT RETAIN