Principal Component Analysis (PCA)

# Principal Components Analysis (PCA): Mental Health Symptom Reduction
# Reducing 15 symptom measures to fewer dimensions

library(psych)        # For PCA, KMO, parallel analysis
library(corrplot)     # Correlation matrix visualization
library(FactoMineR)   # Comprehensive PCA
library(factoextra)   # PCA visualization
library(dplyr)

# Simulate realistic mental health symptom data
set.seed(2025)
n <- 320

# Generate 3 underlying factors
factor1_internalizing <- rnorm(n)  # Depression-anxiety
factor2_stress <- rnorm(n)          # Stress-tension  
factor3_fear <- rnorm(n)            # Fear-panic

# Create 15 observed variables as linear combinations + noise
data <- data.frame(
  # Depression symptoms (load on factor 1)
  dep_mood = 0.85*factor1_internalizing + rnorm(n, 0, 0.3),
  dep_anhedonia = 0.80*factor1_internalizing + rnorm(n, 0, 0.35),
  dep_fatigue = 0.75*factor1_internalizing + 0.3*factor2_stress + rnorm(n, 0, 0.3),
  dep_worthless = 0.82*factor1_internalizing + rnorm(n, 0, 0.32),
  dep_concentration = 0.70*factor1_internalizing + 0.25*factor2_stress + rnorm(n, 0, 0.35),
  
  # Anxiety symptoms (load on factors 1 and 3)
  anx_worry = 0.78*factor1_internalizing + 0.25*factor3_fear + rnorm(n, 0, 0.3),
  anx_restless = 0.72*factor1_internalizing + 0.30*factor2_stress + rnorm(n, 0, 0.35),
  anx_tense = 0.68*factor1_internalizing + 0.35*factor2_stress + rnorm(n, 0, 0.32),
  anx_panic = 0.40*factor1_internalizing + 0.75*factor3_fear + rnorm(n, 0, 0.3),
  anx_fear = 0.35*factor1_internalizing + 0.80*factor3_fear + rnorm(n, 0, 0.28),
  
  # Stress symptoms (load on factor 2)
  stress_overwhelm = 0.30*factor1_internalizing + 0.82*factor2_stress + rnorm(n, 0, 0.3),
  stress_control = 0.25*factor1_internalizing + 0.78*factor2_stress + rnorm(n, 0, 0.32),
  stress_irritable = 0.40*factor1_internalizing + 0.70*factor2_stress + rnorm(n, 0, 0.35),
  stress_pressure = 0.20*factor1_internalizing + 0.80*factor2_stress + rnorm(n, 0, 0.30),
  stress_coping = 0.28*factor1_internalizing + 0.75*factor2_stress + rnorm(n, 0, 0.33)
)

# Standardize to 0-100 scale (simulating normalized scores)
data <- as.data.frame(scale(data) * 15 + 50)
data[data < 0] <- 0
data[data > 100] <- 100

# === STEP 1: Check Assumptions ===

cat("=== Data Summary ===\n")
cat("Sample size:", nrow(data), "\n")
cat("Number of variables:", ncol(data), "\n")
cat("N/p ratio:", round(nrow(data)/ncol(data), 1), "\n\n")

# Check for missing data
cat("Missing data:", sum(is.na(data)), "values\n\n")

# Descriptive statistics
summary(data)

# 1. Check correlations (need adequate correlations)
cor_matrix <- cor(data)
cat("=== Correlation Matrix Summary ===\n")
cat("Range of correlations:", round(min(cor_matrix[cor_matrix!=1]), 2), 
    "to", round(max(cor_matrix[cor_matrix!=1]), 2), "\n\n")

# Visualize correlation matrix
corrplot(cor_matrix, method="color", type="upper", 
         tl.cex=0.7, tl.col="black",
         title="Correlation Matrix: Mental Health Symptoms",
         mar=c(0,0,2,0))

# 2. KMO - Sampling Adequacy
kmo_result <- KMO(data)
cat("=== Kaiser-Meyer-Olkin(KMO) ===\n")
cat("Overall MSA:", round(kmo_result$MSA, 3), "\n")
if (kmo_result$MSA > 0.9) {
  cat("Interpretation: Marvelous - excellent for PCA\n\n")
} else if (kmo_result$MSA > 0.8) {
  cat("Interpretation: Meritorious - very good for PCA\n\n")
} else if (kmo_result$MSA > 0.7) {
  cat("Interpretation: Middling - acceptable for PCA\n\n")
} else if (kmo_result$MSA > 0.6) {
  cat("Interpretation: Mediocre - acceptable but not ideal\n\n")
} else {
  cat("Interpretation: Unacceptable - PCA not recommended\n\n")
}

cat("Individual MSA per variable:\n")
print(round(kmo_result$MSAi, 2))

# 3. Bartlett's Test of Sphericity
bartlett_result <- cortest.bartlett(cor_matrix, n=nrow(data))
cat("\n=== Bartlett's Test of Sphericity ===\n")
cat("Chi-square:", round(bartlett_result$chisq, 2), "\n")
cat("df:", bartlett_result$df, "\n")
cat("p-value:", format(bartlett_result$p.value, scientific=TRUE), "\n")
if (bartlett_result$p.value < 0.05) {
  cat("Interpretation: Significant(p < .05) - correlations adequate for PCA\n\n")
} else {
  cat("Interpretation: Non-significant - variables may be too independent\n\n")
}

# 4. Check for outliers (Mahalanobis distance)
maha_dist <- mahalanobis(data, colMeans(data), cov(data))
outlier_cutoff <- qchisq(0.999, df=ncol(data))
outliers <- which(maha_dist > outlier_cutoff)
cat("Multivariate outliers(p < .001):", length(outliers), "cases\n\n")

# === STEP 2: Determine Number of Components ===

# Method 1: Scree plot
pca_initial <- prcomp(data, center=TRUE, scale.=TRUE)
screeplot(pca_initial, type="lines", 
          main="Scree Plot: Eigenvalues by Component",
          xlab="Component Number", 
          ylab="Eigenvalue(Variance)")
abline(h=1, col="red", lty=2)  # Kaiser criterion line

# Extract eigenvalues
eigenvalues <- pca_initial$sdev^2
cat("=== Eigenvalues(Kaiser Criterion: λ > 1) ===\n")
print(round(eigenvalues, 2))
cat("\nComponents with eigenvalue > 1:", sum(eigenvalues > 1), "\n\n")

# Method 2: Variance explained
var_explained <- eigenvalues / sum(eigenvalues) * 100
cum_var <- cumsum(var_explained)

cat("=== Variance Explained ===\n")
var_table <- data.frame(
  Component = 1:length(eigenvalues),
  Eigenvalue = round(eigenvalues, 2),
  Variance_Pct = round(var_explained, 1),
  Cumulative_Pct = round(cum_var, 1)
)
print(var_table[1:6,])  # First 6 components

# Method 3: Parallel Analysis (BEST METHOD)
parallel <- fa.parallel(data, fa="pc", n.iter=1000, 
                       main="Parallel Analysis: Actual vs Simulated Data")
cat("\n=== Parallel Analysis ===\n")
cat("Suggested number of components:", parallel$ncomp, "\n\n")

# Decision: Retain 3 components based on parallel analysis
n_components <- 3
cat("DECISION: Retaining", n_components, "components\n")
cat("Total variance explained:", round(cum_var[n_components], 1), "%\n\n")

# === STEP 3: Run PCA with Optimal Number of Components ===

pca_result <- principal(data, nfactors=n_components, 
                       rotate="varimax",  # Orthogonal rotation
                       scores=TRUE)

cat("=== PCA Results ===\n")
print(pca_result)

# Component loadings
cat("\n=== Component Loadings(Varimax Rotation) ===\n")
loadings_matrix <- as.data.frame(unclass(pca_result$loadings))
print(round(loadings_matrix, 2))

# Communalities (variance explained per variable)
cat("\n=== Communalities(h²) ===\n")
cat("Proportion of variance in each variable explained by components\n")
print(round(pca_result$communality, 2))
cat("\nVariables with low communality(<0.3):", 
    sum(pca_result$communality < 0.3), "\n\n")

# === STEP 4: Interpret Components ===

cat("=== Component Interpretation ===\n\n")
cat("PC1(Internalizing): Loads high on depression and general anxiety\n")
cat("PC2(Stress): Loads high on stress and tension symptoms\n")
cat("PC3(Fear): Loads high on panic and fear symptoms\n\n")

# === STEP 5: Visualizations ===

# Biplot (variables and observations)
biplot(pca_initial, scale=0, cex=0.6,
       main="PCA Biplot: Variables and Observations")

# Component loadings plot
library(factoextra)
fviz_pca_var(pca_initial, col.var="contrib",
             gradient.cols=c("#00AFBB", "#E7B800", "#FC4E07"),
             repel=TRUE,
             title="PCA Variable Plot: Contributions to PC1-PC2")

# Individual plot (colored by contribution)
fviz_pca_ind(pca_initial, col.ind="cos2",
             gradient.cols=c("#00AFBB", "#E7B800", "#FC4E07"),
             repel=FALSE,
             title="PCA Individual Plot")

# === STEP 6: Save Component Scores ===

component_scores <- as.data.frame(pca_result$scores)
colnames(component_scores) <- c("Internalizing", "Stress", "Fear")

cat("=== Component Scores(First 6 subjects) ===\n")
print(head(component_scores))

# These scores can be used in subsequent analyses
data_with_components <- cbind(data, component_scores)

cat("\n=== Summary ===\n")
cat("Successfully reduced 15 variables to 3 principal components\n")
cat("explaining", round(cum_var[n_components], 1), "% of total variance.\n")
cat("Components represent: Internalizing, Stress, and Fear dimensions.\n")
cat("Component scores saved for use in predictive models.\n")