PCA_significance.R

# Load necessary libraries
library(MASS)      # For generating random matrices
library(RMTstat)   # For Marchenko-Pastur distribution

setwd("/Users/enrique/Dropbox/Ranier_modelos/Archivos_CSV/")

example_data <- read.csv(file="Omission_lateCorregido.csv", header=FALSE, sep=",")

# Function to compute eigenvalues
compute_eigenvalues <- function(data) {
  cov_matrix <- cov(data)
  eigen(cov_matrix)$values
}

# Function to generate Marchenko-Pastur bounds
marchenko_pastur_bounds <- function(p, n) {
  q <- p / n
  lambda_min <- (1 - sqrt(q))^2
  lambda_max <- (1 + sqrt(q))^2
  list(lambda_min = lambda_min, lambda_max = lambda_max)
}

# Main function for PCA significance using Marchenko-Pastur distribution
pca_significance <- function(data, num_permutations = 1000) {
  n <- nrow(data)
  p <- ncol(data)
  
  # Compute actual eigenvalues
  actual_eigenvalues <- compute_eigenvalues(data)
  
  # Generate null eigenvalues through random permutations
  null_eigenvalues <- replicate(num_permutations, {
    permuted_data <- apply(data, 2, sample)
    compute_eigenvalues(permuted_data)
  })
  
  # Flatten null eigenvalues for distribution fitting
  null_eigenvalues <- as.vector(null_eigenvalues)
  
  # Fit the Marchenko-Pastur bounds
  mp_bounds <- marchenko_pastur_bounds(p, n)
  
  # Identify significant eigenvalues
  significant_eigenvalues <- actual_eigenvalues[actual_eigenvalues > mp_bounds$lambda_max]
  
  # Plot results
  hist(null_eigenvalues, breaks = 50, probability = TRUE, 
       col = "lightblue", main = "Eigenvalue Distribution", 
       xlab = "Eigenvalue", ylab = "Density")
  curve(dmp(x, p/n), add = TRUE, col = "red", lwd = 2)
  abline(v = mp_bounds$lambda_min, col = "blue", lty = 2, lwd = 2)
  abline(v = mp_bounds$lambda_max, col = "blue", lty = 2, lwd = 2)
  points(actual_eigenvalues, rep(0, length(actual_eigenvalues)), 
         col = ifelse(actual_eigenvalues > mp_bounds$lambda_max, "red", "black"), 
         pch = 19)
  
  # Return significant eigenvalues and bounds
  list(
    significant_eigenvalues = significant_eigenvalues,
    mp_bounds = mp_bounds,
    actual_eigenvalues = actual_eigenvalues
  )
}

# Example usage
set.seed(123)
#example_data <- matrix(rnorm(1000), nrow = 100, ncol = 10)  # Example data

result <- pca_significance(example_data)

# Print results
cat("Significant Eigenvalues:\n")
print(result$significant_eigenvalues)
cat("Marchenko-Pastur Bounds:\n")
print(result$mp_bounds)