PCA versus Factor Analysis

Purpose

Principal component analysis (PCA) and factor analysis (FA) are often confused, and in practice they are sometimes used interchangeably. Both take a large set of correlated variables and return a few dimensions with loadings and scores. But they rest on different models of where the data come from. This walkthrough:

1. Explains the conceptual and statistical difference between FA and PCA.
2. Explains what psych::principal() is actually doing.
3. Fits a four-dimension solution both ways on the same census data and compares the loadings, the variance explained, and the scores.

It builds on the Factor Analysis walkthrough, which uses the same data and explains how the four-factor solution was chosen.

1. Two different questions

2. Common versus total variance

The practical implication: PCA loadings usually look stronger than FA loadings. PCA counts all of a variable’s variance, including the noise, as something to explain. FA loadings reflect only the variance a variable shares with the others.

3. What psych::principal() is doing

A common snippet looks like this:

pc4 <- principal(d, nfactors = 4, rotate = "varimax", scores = TRUE)
print.psych(pc4, cut = 0.3, sort = TRUE)

This is PCA, even though the argument is named nfactors. psych::principal() extracts principal components, optionally rotates them, and returns component scores. The psych package uses the same argument names for both methods so the output is easy to compare, which is also why the two are easy to confuse.

4. Setup and data

# Install once if needed:
# install.packages(c("tidyverse", "psych", "GPArotation", "pander"))

library(tidyverse)     # data import and wrangling
library(psych)         # principal(), fa()
library(GPArotation)   # rotation methods
library(pander)        # clean tables in the knitted report

# Read local copies if they sit next to this file; otherwise read from the course repo.
repo <- "https://raw.githubusercontent.com/lecy/econometrics/main/lectures2/data/"
data_path <- if (file.exists("census-data.csv"))    "census-data.csv"    else paste0(repo, "census-data.csv")
dd_path   <- if (file.exists("census-data-dd.csv")) "census-data-dd.csv" else paste0(repo, "census-data-dd.csv")

out_dir <- "fa_vs_pca_outputs"
dir.create(out_dir, showWarnings = FALSE)

For comparability, both methods get identical inputs:

1. Median imputation for any missing values (a simple exploratory default)
2. Z-score standardization, so both methods analyze the correlation matrix

dat <- readr::read_csv(data_path, show_col_types = FALSE)
dd  <- readr::read_csv(dd_path,   show_col_types = FALSE)

stopifnot(all(c("LABEL", "VARIABLE") %in% names(dd)))
vars <- dd$LABEL
stopifnot(all(vars %in% names(dat)))

median_impute <- function(x) ifelse(is.na(x), median(x, na.rm = TRUE), x)

X_z <- dat %>%
  dplyr::select(all_of(vars)) %>%
  mutate(across(everything(), ~ median_impute(as.numeric(.x)))) %>%
  scale() %>%
  as.data.frame()

dim(X_z)

## [1] 1000   30

5. Fit PCA and FA with four dimensions

Both use varimax rotation, so each solution has uncorrelated dimensions and is as easy to read as possible.

pca4 <- psych::principal(
  X_z,
  nfactors = 4,           # number of components to keep
  rotate   = "varimax",   # orthogonal rotation
  scores   = TRUE         # return component scores
)

fa4 <- suppressWarnings(psych::fa(
  X_z,
  nfactors = 4,
  fm       = "pa",          # principal axis factoring
  rotate   = "varimax",
  scores   = "regression"   # regression-method factor scores
))

6. Line up the dimensions

Neither method guarantees that its first dimension matches the other method’s first dimension. PCA names its rotated components RC1–RC4 and FA names its factors PA1–PA4, and the numbering doesn’t correspond. Before comparing, match each factor to the component whose scores correlate with it most strongly, and flip the sign if needed. The sign of a loading column is arbitrary in both methods.

cross <- cor(fa4$scores, pca4$scores)
round(cross, 2)

##       RC1  RC4   RC2   RC3
## PA1  0.98 0.04  0.00 -0.02
## PA2 -0.01 0.95 -0.04 -0.01
## PA4  0.00 0.01  0.92  0.00
## PA3  0.00 0.02 -0.02  0.89

# For each FA factor, the PCA component it matches best
match_pc <- colnames(cross)[apply(abs(cross), 1, which.max)]
names(match_pc) <- rownames(cross)
stopifnot(!anyDuplicated(match_pc))
match_sign <- sign(cross[cbind(rownames(cross), match_pc)])
names(match_sign) <- rownames(cross)

# Construct names from the Factor Analysis walkthrough
factor_names <- c(PA1 = "ses", PA2 = "stability", PA4 = "immigrant", PA3 = "aging")
fa_order <- names(factor_names)

tibble(
  construct = factor_names[fa_order],
  fa_factor = fa_order,
  pca_component = match_pc[fa_order],
  score_correlation = round(cross[cbind(fa_order, match_pc[fa_order])], 3)
)

Every factor has a clear partner component. The four dimensions are the same constructs under either method: socioeconomic status, household stability, immigrant communities, and aging population.

7. Compare the loadings side by side

L_fa  <- unclass(fa4$loadings)[, fa_order]
L_pca <- sweep(unclass(pca4$loadings)[, match_pc[fa_order]], 2, match_sign[fa_order], `*`)
colnames(L_fa)  <- paste0("FA_",  factor_names[fa_order])
colnames(L_pca) <- paste0("PCA_", factor_names[fa_order])

loadings_compare <- dd %>%
  dplyr::select(LABEL, VARIABLE) %>%
  left_join(as.data.frame(L_pca) %>% rownames_to_column("LABEL"), by = "LABEL") %>%
  left_join(as.data.frame(L_fa)  %>% rownames_to_column("LABEL"), by = "LABEL")

readr::write_csv(loadings_compare, file.path(out_dir, "loadings_pca_vs_fa.csv"))

Show each construct’s PCA and FA loadings next to each other, for the variables where either method has a loading of at least 0.40 on that construct.

for (k in unname(factor_names[fa_order])) {
  tab <- loadings_compare %>%
    transmute(LABEL, VARIABLE,
              PCA = .data[[paste0("PCA_", k)]],
              FA  = .data[[paste0("FA_",  k)]],
              pca_minus_fa = abs(PCA) - abs(FA)) %>%   # compare sizes, ignoring sign
    filter(pmax(abs(PCA), abs(FA)) >= 0.40) %>%
    arrange(desc(abs(PCA))) %>%
    mutate(across(c(PCA, FA, pca_minus_fa), ~ sprintf("%.2f", .x)))
  cat("\n### ", k, "\n\n", sep = "")
  pander::pandoc.table(tab, justify = "llrrr", split.table = Inf)
}

8. Variance explained and communalities

Both methods report a communality (\(h^2\)) for each variable: the share of its variance captured by the four dimensions. In FA, the rest is uniqueness. In PCA, the rest is simply variance the four components left behind.

comm <- tibble(
  LABEL   = names(fa4$communality),
  h2_pca  = unname(pca4$communality),
  h2_fa   = unname(fa4$communality)
) %>%
  left_join(dd, by = "LABEL") %>%
  relocate(LABEL, VARIABLE)

summarise(comm,
          mean_h2_pca = mean(h2_pca),
          mean_h2_fa  = mean(h2_fa),
          share_pca_higher = mean(h2_pca > h2_fa))

# Variables where FA reports the higher communality
comm %>% filter(h2_fa >= h2_pca) %>% mutate(across(c(h2_pca, h2_fa), ~ round(.x, 2)))

lim <- c(0, 1.05)
plot(comm$h2_fa, comm$h2_pca, xlim = lim, ylim = lim,
     pch = 19, col = "#bd5416", bty = "n",
     xlab = "Communality in factor analysis",
     ylab = "Communality in PCA",
     main = "Communality: PCA versus factor analysis")
abline(0, 1, lty = 2, col = "#6d7681")
text(comm$h2_fa, comm$h2_pca, comm$LABEL, pos = 4, cex = 0.6, col = "#4d565f")

Points above the dashed line are variables where PCA claims more variance than FA, which is 83% of them. The gap tends to be widest for variables with low FA communalities, such as percent 17 and under and female labor force participation. Those are the noisy indicators, where FA attributes most of the variance to uniqueness and PCA keeps it. The few points below the line are the near-duplicates and a handful of very strong indicators.

The same pattern shows up in total variance explained:

round(pca4$Vaccounted[c("SS loadings", "Proportion Var", "Cumulative Var"), ], 2)

##                 RC1  RC4  RC2  RC3
## SS loadings    6.32 4.85 3.56 3.14
## Proportion Var 0.21 0.16 0.12 0.10
## Cumulative Var 0.21 0.37 0.49 0.60

round(fa4$Vaccounted[c("SS loadings", "Proportion Var", "Cumulative Var"), ], 2)

##                 PA1  PA2  PA4  PA3
## SS loadings    6.13 4.44 3.27 2.79
## Proportion Var 0.20 0.15 0.11 0.09
## Cumulative Var 0.20 0.35 0.46 0.55

9. Compare the scores

Both methods produce a score on each dimension for every neighborhood. How similar are they, and how independent is each set?

scores_fa <- as.data.frame(fa4$scores)[, fa_order]
names(scores_fa) <- paste0("FA_", factor_names[fa_order])

scores_pca <- sweep(as.data.frame(pca4$scores)[, match_pc[fa_order]], 2, match_sign[fa_order], `*`) %>%
  as.data.frame()
names(scores_pca) <- paste0("PCA_", factor_names[fa_order])

# Within-method correlations
round(cor(scores_pca), 3)

##               PCA_ses PCA_stability PCA_immigrant PCA_aging
## PCA_ses             1             0             0         0
## PCA_stability       0             1             0         0
## PCA_immigrant       0             0             1         0
## PCA_aging           0             0             0         1

round(cor(scores_fa), 3)

##              FA_ses FA_stability FA_immigrant FA_aging
## FA_ses        1.000       -0.002        0.000    0.013
## FA_stability -0.002        1.000        0.026   -0.031
## FA_immigrant  0.000        0.026        1.000    0.007
## FA_aging      0.013       -0.031        0.007    1.000

# Matched pairs across methods
round(diag(cor(scores_pca, scores_fa)), 3)

## [1] 0.983 0.945 0.922 0.886

readr::write_csv(scores_pca, file.path(out_dir, "scores_pca4.csv"))
readr::write_csv(scores_fa,  file.path(out_dir, "scores_fa4.csv"))
readr::write_csv(
  as.data.frame(round(cor(scores_pca, scores_fa), 3)) %>% rownames_to_column("PCA_component"),
  file.path(out_dir, "score_correlations_pca_vs_fa.csv")
)

Three things to notice:

• PCA component scores are exactly uncorrelated. They are exact projections onto orthogonal axes.
• FA factor scores are only approximately uncorrelated. The factors are orthogonal, but the scores are estimates of them, and estimation error introduces small correlations.
• The matched scores are highly correlated across methods (roughly 0.9 or higher). For a large, well-structured set of variables, the two methods agree on which dimensions exist. They disagree about how much of each variable is signal.

10. When to use which

Outputs

Knitting this file writes four files to fa_vs_pca_outputs/:

• loadings_pca_vs_fa.csv: aligned PCA and FA loadings for every variable
• scores_pca4.csv: aligned PCA component scores
• scores_fa4.csv: FA factor scores
• score_correlations_pca_vs_fa.csv: correlations between the two sets of scores