##########################################################
# Factor Analysis in R - Chestnut Ridge Example          #
# Retailer Brand Audit from Survey Data                  #
##########################################################

# This script reproduces the full Chestnut Ridge retailer brand audit example in R:
#   Step 1. Look at the average survey responses by retailer
#   Step 2. Determine the number of factors (eigenvalues and scree plot)
#   Step 3. Run a factor analysis with four factors (default varimax rotation)
#   Step 4. Rerun the factor analysis with an oblique rotation (heat map of loadings)
#   Step 5. Name the factors and score each retailer on the factors
#   Step 6. Heat map of factor scores by retailer
# Run it from top to bottom. Each step prints its result to the console or Plots pane.

## Install Packages (if needed - only installs packages that are missing)
packages <- c("dplyr", "tidyr", "ggplot2", "GPArotation")
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) install.packages(new_packages)

## Load Packages
library(dplyr)        # data manipulation
library(tidyr)        # reshaping tables for the heat maps
library(ggplot2)      # plots and heat maps
library(GPArotation)  # oblique (oblimin) factor rotation

## Set Seed
# Setting the seed ensures your results are the same as the results in this example.
set.seed(1)

# Import Data
retailer_survey <- read.csv(file.choose()) ## Choose retailer_survey.csv file

# Take a first look at the data
# 600 rows: 100 respondents rating each of 6 retailers (Chestnut Ridge = CR, and A-E)
# on 12 statements, each measured on a 1-7 scale
head(retailer_survey)
table(retailer_survey$Retailer)


#####################################################
# Step 1: Average Survey Responses by Retailer      #
#####################################################

# Average response to each of the 12 statements for each retailer
retailer_means <- retailer_survey %>%
  group_by(Retailer) %>%
  summarise(across(everything(), mean))

cat("\nAverage Survey Responses by Retailer\n")
print(as.data.frame(mutate(retailer_means, across(-Retailer, ~ round(.x, 2)))), row.names = FALSE)


#####################################################
# Step 2: Determine the Number of Factors           #
#####################################################

# Factor analysis uses only the 12 statements, so remove the retailer column
retailer_factors <- select(retailer_survey, -Retailer)

# Eigenvalues of the correlation matrix
eigenvalues <- eigen(cor(retailer_factors))$values
eigenvalues

# Latent root (eigenvalue) criterion: keep factors with an eigenvalue of at least 1
sum(eigenvalues >= 1)

# Percentage of variance criterion: share of the variance explained by each factor
variance_explained <- data.frame(factor = 1:length(eigenvalues),
                                 eigenvalue = eigenvalues,
                                 pct_variance = eigenvalues / sum(eigenvalues),
                                 cumulative_pct = cumsum(eigenvalues) / sum(eigenvalues))
print(mutate(variance_explained, across(-factor, ~ round(.x, 3))), row.names = FALSE)

# Scree plot of the eigenvalues
ggplot(variance_explained, aes(x = factor, y = eigenvalue)) +
  geom_line(color = "blue") +
  geom_point(size = 2) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
  scale_x_continuous(breaks = 1:length(eigenvalues)) +
  labs(title = "Scree Plot of Eigenvalues",
       subtitle = "Dashed line marks an eigenvalue of 1",
       x = "Factor", y = "Eigenvalue") +
  theme_minimal()

# Four eigenvalues are greater than 1 (3.27, 2.90, 2.67, and 2.05), and together
# they explain about 91% of the variance, so we use a four-factor solution.


#####################################################
# Step 3: Factor Analysis with Four Factors         #
#####################################################

# Run factor analysis with 4 factors (default varimax rotation)
# Loadings close to 0 (below 0.1) are left blank in the printed output
factanal(retailer_factors, factors = 4)


#####################################################
# Step 4: Factor Analysis with Oblique Rotation     #
#####################################################

# Rerun the factor analysis with an oblique rotation and save the factor scores
retailer_fa <- factanal(retailer_factors, factors = 4,
                        rotation = "oblimin", scores = "Bartlett")

# Factor loadings (all loadings are shown, including those close to 0)
print(retailer_fa$loadings, cutoff = 0, digits = 2)

# Put the loadings in a table for plotting
loadings_table <- as.data.frame(unclass(retailer_fa$loadings)) %>%
  mutate(statement = rownames(.)) %>%
  pivot_longer(-statement, names_to = "factor", values_to = "loading")

# List the statements so that statements on the same factor are next to each other
statement_order <- c("Latest", "Trends", "Stylish", "Quality", "Last", "Fit",
                     "Satisfied", "Purchase", "Recommend", "Value", "Bargain", "Worth")

# Heat map of factor loadings (darker red means a stronger factor loading)
ggplot(loadings_table, aes(x = factor, y = factor(statement, levels = rev(statement_order)),
                           fill = loading)) +
  geom_tile(color = "white") +
  geom_text(aes(label = sprintf("%.2f", loading)), size = 3.5) +
  scale_fill_distiller(palette = "Reds", direction = 1, name = "Loading") +
  labs(title = "Factor Loadings from Survey", x = NULL, y = NULL) +
  theme_minimal() +
  theme(panel.grid = element_blank())


#####################################################
# Step 5: Name the Factors and Score the Retailers  #
#####################################################

# The three statements with the highest loadings on each factor (Table 11.4)
loadings_table %>%
  group_by(factor) %>%
  slice_max(loading, n = 3) %>%
  summarise(highest_loadings = paste(statement, collapse = ", "))

# Name each factor to capture the essence of the statements that load on it
factor_names <- c("Innovative", "High Quality", "Loyalty", "Good Value")

# Factor scores for each response, with the retailer that was rated
retailer_scores <- data.frame(retailer_fa$scores)
names(retailer_scores) <- factor_names
retailer_scores$Retailer <- retailer_survey$Retailer

# Average factor score for each retailer
# Factor scores are standardized (mean 0, standard deviation 1), so a positive score
# means the retailer is perceived above average on that factor
retailer_fa_mean <- retailer_scores %>%
  group_by(Retailer) %>%
  summarise(across(all_of(factor_names), mean))

cat("\nAverage Factor Scores by Retailer\n")
print(as.data.frame(mutate(retailer_fa_mean, across(-Retailer, ~ round(.x, 3)))), row.names = FALSE)


#####################################################
# Step 6: Heat Map of Factor Scores by Retailer     #
#####################################################

# Put the average factor scores in a table for plotting
scores_table <- retailer_fa_mean %>%
  pivot_longer(-Retailer, names_to = "factor", values_to = "score")

# Heat map of factor scores by retailer (darker blue means a higher factor score)
ggplot(scores_table, aes(x = factor(factor, levels = factor_names),
                         y = factor(Retailer, levels = rev(sort(unique(Retailer)))),
                         fill = score)) +
  geom_tile(color = "white") +
  geom_text(aes(label = sprintf("%.2f", score)), size = 3.5) +
  scale_fill_distiller(palette = "Blues", direction = 1, name = "Factor Score") +
  labs(title = "Factor Score by Retailer", x = NULL, y = "Retailer") +
  theme_minimal() +
  theme(panel.grid = element_blank())