Реализация параллельного анализа Хорна для PCA на ковариационных матрицах

# Загрузка необходимых библиотек
library(nFactors)
library(psych)

# Импорт ваших данных
data <- read.xlsx("data.xlsx")

# Убедитесь, что ваши данные правильно масштабированы
# Важно: НЕ стандартизируйте данные для PCA ковариационной матрицы
n_subjects <- nrow(data)
n_variables <- ncol(data)

# Вычисление ковариационной матрицы
cov_mat <- cov(data, use = "pairwise.complete.obs")

# Вычисление собственных значений из ваших фактических данных
eigen_observed <- eigen(cov_mat)$values

# Выполните параллельный анализ для ковариационной матрицы
results <- parallel(
  var = n_variables,
  subject = n_subjects,
  rep = 1000,          # Количество репликаций
  quantile = 0.95,     # Порог 95-го процентиля
  model = "components", # Модель PCA
  sd = sqrt(diag(cov_mat)) # Критически важно: использовать фактические стандартные отклонения
)

# Извлечение смоделированных собственных значений
simulated_eigen_mean <- results$eigen$mevpea
simulated_eigen_quantile <- results$eigen$qevpea

# Создание фрейма данных для сравнения
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
)

# Просмотр результатов
print(comparison)

library(paran)

# Выполнение параллельного анализа
pa_result <- paran(
  cov_mat,
  n.iter = 1000,
  cent = 0.05,
  quantile = TRUE
)

# Построение графика результатов
plot(pa_result)

# Получение рекомендуемого количества компонент
pa_result$nfact

# Ручная реализация параллельного анализа
set.seed(123)  # Для воспроизводимости
n_sim <- 1000
sim_eigen <- matrix(NA, nrow = n_variables, ncol = n_sim)

for (i in 1:n_sim) {
  # Генерация случайных данных с той же ковариационной структурой
  random_data <- matrix(rnorm(n_subjects * n_variables), 
                       nrow = n_subjects, ncol = n_variables)
  
  # Масштабирование для соответствия исходным дисперсиям (критически важно для ковариационной матрицы)
  random_data_scaled <- sweep(random_data, 2, 
                             sqrt(diag(cov_mat)), `*`)
  
  # Вычисление ковариационной матрицы и собственных значений
  cov_sim <- cov(random_data_scaled)
  sim_eigen[, i] <- eigen(cov_sim)$values
}

# Вычисление среднего и квантилей
sim_mean <- apply(sim_eigen, 1, mean)
sim_95th <- apply(sim_eigen, 1, quantile, probs = 0.95)

# Сравнение с наблюдаемыми собственными значениями
comparison <- data.frame(
  Component = 1:n_variables,
  Observed = eigen_observed,
  Simulated_Mean = sim_mean,
  Simulated_95th = sim_95th,
  Retain = eigen_observed > sim_95th
)

# Полный параллельный анализ Хорна для PCA ковариационной матрицы

# Загрузка библиотек
library(nFactors)
library(psych)
library(openxlsx)

# Установка seed для воспроизводимости
set.seed(123)

# Импорт данных
data <- read.xlsx("data.xlsx")

# Базовая информация о данных
cat("Количество наблюдений:", nrow(data), "\n")
cat("Количество переменных:", ncol(data), "\n")

# Вычисление ковариационной матрицы (НЕ стандартизируйте для PCA ковариации)
cov_mat <- cov(data, use = "pairwise.complete.obs")

# Вычисление наблюдаемых собственных значений
eigen_observed <- eigen(cov_mat)$values

# Выполнение параллельного анализа
results <- parallel(
  var = ncol(data),
  subject = nrow(data),
  rep = 2000,          # Увеличенное количество репликаций для стабильности
  quantile = 0.95,     # 95-й процентиль
  model = "components",
  sd = sqrt(diag(cov_mat))  # Критический параметр для ковариационных матриц
)

# Извлечение смоделированных собственных значений
sim_mean <- results$eigen$mevpea
sim_95th <- results$eigen$qevpea

# Создание фрейма данных для сравнения
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
)

# Отображение результатов
print("Сравнение собственных значений:")
print(comparison)

# Определение количества компонент для сохранения
n_components <- sum(eigen_observed > sim_95th)
cat("\nРекомендуемое количество компонент для сохранения:", n_components, "\n")

# Создание визуализации
plot(1:length(eigen_observed), eigen_observed, 
     type = "b", pch = 19, col = "blue",
     xlab = "Номер компоненты", 
     ylab = "Собственное значение",
     main = "Параллельный анализ Хорна - Ковариационная матрица")
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("Наблюдаемые", "Средние смоделированные", "95-й процентиль смоделированных"),
       col = c("blue", "red", "green"),
       pch = c(19, 17, 15))

# Дополнительные диагностические графики
par(mfrow = c(2, 2))
plotuScree(eigen_observed, main = "График крутизны")
plotuScree(sim_95th, main = "95-й процентиль смоделированных")
hist(eigen_observed, main = "Распределение наблюдаемых собственных значений")
hist(sim_95th, main = "Распределение 95-го процентиля смоделированных")
par(mfrow = c(1, 1))

Компонента 1: Наблюдаемое = 5.2, С��елированное_95-е = 1.8 → СОХРАНИТЬ
Компонента 2: Наблюдаемое = 2.1, С��елированное_95-е = 1.5 → СОХРАНИТЬ
Компонента 3: Наблюдаемое = 1.2, С��елированное_95-е = 1.3 → НЕ СОХРАНЯТЬ