Factor Analysis

Purpose

The data contain 30 census variables for 1,000 neighborhoods (census tracts). Each one plausibly describes some part of neighborhood health, stability, or robustness, and many of them are highly correlated. Median household income, per capita income, and the share of adults with a college degree are not really three separate facts about a neighborhood.

Factor analysis asks whether a smaller number of latent constructs could be generating the correlations we see. The goals of this walkthrough are to:

Check that the variables share enough correlation to be worth factoring.
Choose the number of factors with a scree plot and parallel analysis.
Interpret the factor loadings and give each factor a name.
Extract factor scores, the new composite measures, for use in later models or maps.

What factor analysis assumes

Every observed variable is a weighted combination of a few shared latent factors, plus a part that is unique to that variable:

\[ X = \Lambda F + \epsilon \]

$X$: the observed variables (the 30 census measures, standardized)
$F$: the latent factors (the constructs we want to recover)
$\Lambda$: the loadings, how strongly each variable reflects each factor
$\epsilon$: the uniqueness, variance in each variable that the factors don’t explain, including measurement error

This is the multi-item version of the measurement model $M = T + e$ from the lecture notes. The factor scores it produces are weighted composites that give more weight to the variables that measure a construct most cleanly. See Weighted Composite Indices for why that beats a simple average.

Three outputs matter most:

Loadings: which variables define each factor. This is how you interpret and name it.
Communality ($h^2$): the share of a variable’s variance explained by all the factors.
Factor scores: each neighborhood’s estimated value on each factor.

1. Setup

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

library(tidyverse)     # data import and wrangling
library(psych)         # fa(), fa.parallel(), KMO(), cortest.bartlett(), alpha()
library(GPArotation)   # rotation methods used by psych::fa()
library(corrplot)      # correlation heat map
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 <- "factor_analysis_outputs"
dir.create(out_dir, showWarnings = FALSE)

psych does the heavy lifting:

fa.parallel() draws the scree plot and runs parallel analysis
fa() fits the factor model
KMO() and cortest.bartlett() check whether factoring is appropriate at all

2. Read the data and the dictionary

The data dictionary census-data-dd.csv has two columns: LABEL, the variable name in census-data.csv, and VARIABLE, a readable description.

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

missing_vars <- setdiff(vars, names(dat))
if (length(missing_vars) > 0) {
  stop("These dictionary variables are missing from census-data.csv: ",
       paste(missing_vars, collapse = ", "))
}

X <- dat %>%
  dplyr::select(all_of(vars)) %>%
  mutate(across(everything(), as.numeric))

dim(X)

## [1] 1000   30

pander(dd, justify = "left")

LABEL	VARIABLE
pnhwht12	Percent white, non-Hispanic
pnhblk12	Percent black, non-Hispanic
phisp12	Percent Hispanic
pntv12	Percent Native American race
pfb12	Percent foreign born
polang12	Percent speaking other language at home, age 5 plus
phs12	Percent with high school degree or less
pcol12	Percent with 4-year college degree or more
punemp12	Percent unemployed
pflabf12	Percent female labor force participation
pprof12	Percent professional employees
pmanuf12	Percent manufacturing employees
pvet12	Percent veteran
psemp12	Percent self-employed
hinc12	Median HH income, total
incpc12	Per capita income
ppov12	Percent in poverty, total
pown12	Percent owner-occupied units
pvac12	Percent vacant units
pmulti12	Percent multi-family units
mrent12	Median rent
mhmval12	Median home value
p30old12	Percent structures more than 30 years old
p10yrs12	Percent HH in neighborhood 10 years or less
p18und12	Percent 17 and under, total
p60up12	Percent 60 and older, total
p75up12	Percent 75 and older, total
pmar12	Percent currently married, not separated
pwds12	Percent widowed, divorced and separated
pfhh12	Percent female-headed families with children

The correlation structure

Before fitting anything, look at the correlations. order = "hclust" reorders the variables so that clusters of highly correlated measures sit next to each other. Those blocks are the factors we expect to find.

corr_matrix <- cor(X, use = "complete.obs")

corrplot(
  corr_matrix,
  order = "hclust",    # group variables with similar correlation profiles
  tl.col = "black",    # label color
  tl.cex = 0.75        # label size
)

3. Standardize the variables

Why standardize?

Factor analysis is fit to the correlation matrix, which is the covariance matrix of z-scored variables. Standardizing puts median home value (in dollars) and percent unemployed (in percentage points) on the same scale, so each variable contributes through its relationships rather than its units.

Here are four raw variables:

dsub <- X[c("pown12", "pvet12", "pnhwht12", "polang12")]
pairs(dsub, lower.panel = panel.smooth, upper.panel = panel.cor)

Missing values

These data happen to be complete, but the code below is a safeguard for other versions of the file. It uses median imputation for each variable, a simple and transparent default for exploratory work. For published inference, consider multiple imputation or full information maximum likelihood instead.

sum(is.na(X))   # number of missing values

## [1] 0

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

X_imp <- X %>% mutate(across(everything(), median_impute))

# scale() subtracts each column's mean and divides by its standard deviation
X_z <- as.data.frame(scale(X_imp))

The same four variables after standardizing:

dsub2 <- X_z[c("pown12", "pvet12", "pnhwht12", "polang12")]
histx(dsub2)

pairs(dsub2, lower.panel = panel.smooth, upper.panel = panel.cor)

Compare the two scatterplot matrices. The axes changed, but every correlation is identical. Z-scoring is a linear rescaling, so it cannot change a correlation. That is exactly why factor analysis on standardized data and factor analysis on the correlation matrix are the same thing.

A sensitivity check for outliers

Several of these variables are heavily skewed. When large outliers are present, analysts sometimes winsorize: pull every value below the 5th percentile up to the 5th percentile, and every value above the 95th down to the 95th. This limits the leverage of extreme tracts. It is usually run as a sensitivity analysis, to confirm the results are not driven by a few outliers, rather than as the main analysis.

winsorize_z <- function(x, lower = 0.05, upper = 0.95) {
  x[is.na(x)] <- median(x, na.rm = TRUE)
  lo <- quantile(x, lower)
  hi <- quantile(x, upper)
  x[x < lo] <- lo
  x[x > hi] <- hi
  as.numeric(scale(x))
}

dsub3 <- as.data.frame(lapply(dsub, winsorize_z))
histx(dsub3)

pairs(dsub3, lower.panel = panel.smooth, upper.panel = panel.cor)

Winsorizing does change the correlations, because it is not a linear transformation. The rest of this walkthrough uses the ordinary z-scores in X_z.

4. Is the data factorable?

Two standard checks:

Kaiser-Meyer-Olkin (KMO) measures how much of the correlation between each pair of variables is shared with the other variables. Values above 0.70 are generally good enough for exploratory factor analysis. KMO is reported overall and for each variable.
Bartlett’s test of sphericity tests whether the correlation matrix is different from an identity matrix, one where every variable is uncorrelated with every other. A small p-value means there is shared correlation to factor.

R <- cor(X_z)

kmo  <- psych::KMO(R)
bart <- psych::cortest.bartlett(R, n = nrow(X_z))

round(kmo$MSA, 3)             # overall KMO

## [1] 0.785

sort(round(kmo$MSAi, 2))[1:6] # the six lowest item-level values

##   pntv12 p18und12 pnhblk12  p60up12 p10yrs12 pflabf12 
##     0.11     0.49     0.60     0.64     0.65     0.66

bart

## $chisq
## [1] 29379.03
## 
## $p.value
## [1] 0
## 
## $df
## [1] 435

The overall KMO of 0.78 and Bartlett’s p-value near zero both say the variables are suitable for factoring.

Watch the item-level values. Percent Native American (pntv12) has a KMO of just 0.11. Variables below 0.50 share very little with the rest of the set and are candidates for removal. We keep all 30 here so the example matches the measurement lab, but we will see this variable load on nothing.

5. How many factors?

Scree plot

A scree plot shows the eigenvalues in descending order. Each eigenvalue is the amount of variance a factor accounts for. Look for the elbow, the point where adding another factor buys much less than the one before.

Parallel analysis

Parallel analysis makes the elbow less subjective. It simulates random data with the same number of rows and columns, computes the eigenvalues of that random data, and keeps the factors whose observed eigenvalues beat the random benchmark.

psych::fa.parallel() arguments:

X_z: the standardized data
fm = "pa": principal axis factoring, a robust default for exploratory factor analysis
fa = "fa": report factor analysis eigenvalues only, not principal components
n.iter = 100: number of random data sets to simulate

set.seed(1234)  # parallel analysis uses random draws

fa_par <- psych::fa.parallel(
  X_z,
  fm = "pa",
  fa = "fa",
  n.iter = 100,
  show.legend = TRUE,
  main = "Scree plot and parallel analysis"
)

## Parallel analysis suggests that the number of factors =  8  and the number of components =  NA

fa_par$nfact                 # number of factors suggested by parallel analysis

## [1] 8

round(fa_par$fa.values[1:8], 2)

## [1] 7.44 4.13 1.93 1.35 1.06 0.59 0.38 0.26

Choose the number of factors

Parallel analysis suggests 8 factors. With 1,000 observations, even small factors beat random noise, so parallel analysis tends to recommend more factors as samples grow. The scree plot tells a different story. The first four eigenvalues (7.4, 4.1, 1.9, 1.3) are large, and after the fifth (1.1) the curve is essentially flat.

The number of factors is ultimately a judgment call that balances statistical evidence against interpretability. Here we use four factors. Each is large and each has a clear substantive meaning, and together they account for more than half of the variance in the 30 variables. Try n_factors <- 5 or 6 and see whether the extra factors are worth naming.

n_factors <- 4

6. Fit the factor model

psych::fa() fits the model. The key arguments:

nfactors = n_factors: the number of factors chosen above
fm = "pa": principal axis factoring, which makes fewer distributional assumptions than maximum likelihood
rotate = "varimax": an orthogonal rotation, which keeps the factors uncorrelated so the scores work as independent features
scores = "regression": estimate factor scores with the regression method

Rotation doesn’t change how well the model fits. It spins the factors so that each variable loads strongly on as few factors as possible, which makes the solution easier to read. If you believe the constructs should be correlated, which is plausible for neighborhoods, use an oblique rotation such as rotate = "oblimin" instead.

fa_fit <- psych::fa(
  X_z,
  nfactors = n_factors,
  fm       = "pa",
  rotate   = "varimax",
  scores   = "regression"
)

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected.  Examine the results carefully

About that warning. An ultra-Heywood case means at least one variable has an estimated communality at or slightly above 1, which is impossible for a true variance share. Two variables hit that ceiling here (look for h2 values of 1.00 and negative u2 values in the printout below). Percent speaking another language at home (polang12) is correlated 0.92 with percent foreign born, and percent 60 and older (p60up12) is correlated 0.87 with percent 75 and older. When two indicators are nearly duplicates, a factor can be defined almost entirely by them. The solution is still usable for exploration, and Section 8 checks that it doesn’t depend on the estimation method. For a final instrument, drop or combine redundant indicators.

The setup chunk hides warnings to keep the report tidy, which would have hidden this one. This chunk turns warnings back on with warning=TRUE, because the warnings from a model fit should always be read.

print(fa_fit, digits = 2, cut = 0.30, sort = TRUE)

## Factor Analysis using method =  pa
## Call: psych::fa(r = X_z, nfactors = n_factors, rotate = "varimax", 
##     scores = "regression", fm = "pa")
## Standardized loadings (pattern matrix) based upon correlation matrix
##          item   PA1   PA2   PA4   PA3     h2       u2 com
## pcol12      8  0.94                   0.9119  8.8e-02 1.1
## incpc12    16  0.86                   0.8129  1.9e-01 1.2
## pprof12    11  0.84                   0.7544  2.5e-01 1.1
## phs12       7 -0.82                   0.7747  2.3e-01 1.3
## mhmval12   22  0.76                   0.6521  3.5e-01 1.2
## hinc12     15  0.75  0.45             0.8280  1.7e-01 1.9
## mrent12    21  0.66        0.31       0.5473  4.5e-01 1.5
## punemp12    9 -0.40 -0.33             0.2889  7.1e-01 2.2
## psemp12    14  0.40  0.37             0.3327  6.7e-01 2.4
## pntv12      4                         0.0076  9.9e-01 1.6
## pmar12     28  0.32  0.88             0.8787  1.2e-01 1.3
## pown12     18        0.84             0.7533  2.5e-01 1.1
## pmulti12   20       -0.74             0.6612  3.4e-01 1.5
## pnhwht12    1        0.56 -0.53       0.6566  3.4e-01 2.5
## pfhh12     30 -0.41 -0.55             0.4861  5.1e-01 1.9
## ppov12     17 -0.51 -0.55             0.5738  4.3e-01 2.1
## pnhblk12    2       -0.52             0.3567  6.4e-01 1.6
## pvet12     13        0.37 -0.35  0.32 0.3663  6.3e-01 3.0
## p10yrs12   24       -0.31             0.1689  8.3e-01 2.6
## pmanuf12   12                         0.1545  8.5e-01 2.4
## p30old12   23                         0.1183  8.8e-01 2.8
## polang12    6              0.98       1.0031 -3.1e-03 1.1
## pfb12       5              0.86       0.8495  1.5e-01 1.3
## phisp12     3              0.84       0.7415  2.6e-01 1.1
## p60up12    26        0.33        0.92 1.0001 -5.7e-05 1.3
## p75up12    27                    0.85 0.7616  2.4e-01 1.1
## pwds12     29 -0.35              0.53 0.4424  5.6e-01 2.1
## pflabf12   10                   -0.42 0.2593  7.4e-01 1.8
## p18und12   25                   -0.41 0.3329  6.7e-01 2.9
## pvac12     19                    0.35 0.1571  8.4e-01 1.6
## 
##                        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
## Proportion Explained  0.37 0.27 0.20 0.17
## Cumulative Proportion 0.37 0.64 0.83 1.00
## 
## Mean item complexity =  1.8
## Test of the hypothesis that 4 factors are sufficient.
## 
## df null model =  435  with the objective function =  29.73 with Chi Square =  29379.03
## df of  the model are 321  and the objective function was  9.34 
## 
## The root mean square of the residuals (RMSR) is  0.07 
## The df corrected root mean square of the residuals is  0.08 
## 
## The harmonic n.obs is  1000 with the empirical chi square  1944.95  with prob <  5.4e-230 
## The total n.obs was  1000  with Likelihood Chi Square =  9203.59  with prob <  0 
## 
## Tucker Lewis Index of factoring reliability =  0.583
## RMSEA index =  0.166  and the 90 % confidence intervals are  0.164 0.169
## BIC =  6986.2
## Fit based upon off diagonal values = 0.95
## Measures of factor score adequacy             
##                                                    PA1  PA2 PA4  PA3
## Correlation of (regression) scores with factors   0.98 0.97   1 1.00
## Multiple R square of scores with factors          0.96 0.93   1 0.99
## Minimum correlation of possible factor scores     0.93 0.87   1 0.99

The bottom of the printout reports fit statistics. The Tucker-Lewis index of 0.58 and RMSEA of 0.17 show that four factors leave plenty of correlation unexplained. That is expected when 30 heterogeneous census variables are compressed into four dimensions. The goal here is interpretable constructs, not a close-fitting confirmatory model.

7. Interpret the loadings

Loadings describe how strongly each variable relates to each factor. Rules of thumb:

|loading| ≥ 0.30: potentially meaningful
|loading| ≥ 0.40: strong, in social science contexts
A variable with large loadings on more than one factor is cross-loading, which makes interpretation harder.

The sign matters. A negative loading means the variable moves in the opposite direction from the factor.

print(fa_fit$loadings, cutoff = 0.30, sort = TRUE)

## 
## Loadings:
##          PA1    PA2    PA4    PA3   
## phs12    -0.825                     
## pcol12    0.940                     
## pprof12   0.841                     
## hinc12    0.752  0.447              
## incpc12   0.859                     
## mrent12   0.659         0.314       
## mhmval12  0.764                     
## pnhwht12         0.556 -0.525       
## pnhblk12        -0.518              
## ppov12   -0.510 -0.548              
## pown12           0.844              
## pmulti12        -0.736              
## pmar12    0.317  0.878              
## pfhh12   -0.411 -0.552              
## phisp12                 0.835       
## pfb12                   0.857       
## polang12                0.977       
## p60up12          0.327         0.924
## p75up12                        0.850
## pwds12   -0.348                0.526
## pntv12                              
## punemp12 -0.401 -0.326              
## pflabf12                      -0.423
## pmanuf12                            
## pvet12           0.373 -0.346  0.324
## psemp12   0.396  0.374              
## pvac12                         0.348
## p30old12                            
## p10yrs12        -0.307              
## p18und12                      -0.414
## 
##                  PA1   PA2   PA4   PA3
## SS loadings    6.134 4.439 3.271 2.789
## Proportion Var 0.204 0.148 0.109 0.093
## Cumulative Var 0.204 0.352 0.461 0.554

Name the factors

Reading down each column gives the factors their names:

PA1, Socioeconomic status: college degrees, professional jobs, income, rents, and home values load high; the share with a high school degree or less loads strongly negative.
PA2, Household stability: homeownership and marriage load high; multi-family housing, female-headed families, and poverty load negative.
PA4, Immigrant communities: foreign born, speaking another language at home, and Hispanic ethnicity.
PA3, Aging population: residents 60 and older, 75 and older, and widowed or divorced load high; children and female labor force participation load negative.

Names are claims, not facts. Percent white loads positively on the stability factor and percent Black loads negatively, and percent white also loads negatively on the immigrant factor. A factor labeled “stability” is partly tracking racial composition. Whether that reflects residential segregation, differences in access to homeownership, or something else is a question about validity that the loadings alone cannot answer.

psych numbers factors in the order they were extracted, and after rotation it sorts the columns by variance explained. That is why PA4 appears before PA3.

factor_names <- c(
  PA1 = "ses",
  PA2 = "stability",
  PA4 = "immigrant",
  PA3 = "aging"
)
stopifnot(setequal(names(factor_names), colnames(fa_fit$loadings)))

See the structure

The printed table is precise but hard to scan. A loadings heatmap shows the same numbers as a picture. Each row is a variable and each column a factor. Blue cells load positively, orange cells load negatively, and pale cells are near zero. Rows are grouped by the factor each variable loads on most strongly, the same order as sort = TRUE above, except that variables with no loading of at least 0.30 are collected at the bottom. The bar chart on the right shows each variable’s communality, how much of it the four factors explain.

Three things jump out that are easy to miss in the printout:

Clean blocks down the diagonal are the four constructs.
Cross-loadings show up as a second colored cell in the same row, such as median household income on both socioeconomic status and stability.
Short communality bars flag variables the factors barely explain, like percent Native American and older housing stock.

library(patchwork)   # combine two ggplots side by side

L_long <- unclass(fa_fit$loadings)[, names(factor_names)] %>%
  as.data.frame() %>%
  setNames(factor_names) %>%
  rownames_to_column("LABEL") %>%
  left_join(dd, by = "LABEL") %>%
  mutate(
    communality = fa_fit$communality[LABEL],
    # dominant factor = column with the largest absolute loading
    primary = factor_names[apply(abs(unclass(fa_fit$loadings)[LABEL, names(factor_names)]), 1, which.max)],
    peak = apply(abs(unclass(fa_fit$loadings)[LABEL, names(factor_names)]), 1, max),
    # variables with no loading of at least 0.30 get their own group at the bottom
    primary = ifelse(peak < 0.30, "none", primary),
    primary = factor(primary, levels = c(factor_names, "none"))
  ) %>%
  arrange(primary, desc(peak)) %>%
  mutate(VARIABLE = factor(VARIABLE, levels = rev(unique(VARIABLE))))

heat_data <- L_long %>%
  pivot_longer(all_of(unname(factor_names)), names_to = "factor", values_to = "loading") %>%
  mutate(factor = factor(factor, levels = factor_names))

# horizontal rules between the groups of variables
breaks_y <- L_long %>% count(primary) %>% mutate(y = nrow(L_long) - cumsum(n) + 0.5) %>%
  slice(-n())

factor_labels <- c(ses = "Socioeconomic\nstatus", stability = "Household\nstability",
                   immigrant = "Immigrant\ncommunities", aging = "Aging\npopulation")

p_heat <- ggplot(heat_data, aes(x = factor, y = VARIABLE, fill = loading)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(aes(label = ifelse(abs(loading) >= 0.30, sprintf("%.2f", loading), ""),
                color = abs(loading) > 0.60),
            size = 3.2, show.legend = FALSE) +
  geom_hline(yintercept = breaks_y$y, color = "#4d565f", linewidth = 0.4) +
  scale_fill_gradient2(low = "#bd5416", mid = "#f0efec", high = "#245b8a",
                       midpoint = 0, limits = c(-1, 1),
                       name = "Loading", breaks = c(-1, -0.5, 0, 0.5, 1)) +
  scale_color_manual(values = c(`FALSE` = "#1f2429", `TRUE` = "white")) +
  scale_x_discrete(position = "top", labels = factor_labels) +
  labs(x = NULL, y = NULL) +
  theme_minimal(base_size = 12) +
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(face = "bold", color = "#1f2429"),
        axis.text.y = element_text(color = "#1f2429"),
        legend.position = "bottom",
        legend.key.width = unit(1.6, "cm"))

p_comm <- ggplot(L_long, aes(x = pmin(communality, 1), y = VARIABLE)) +
  geom_col(fill = "#6d7681", width = 0.7) +
  scale_x_continuous(limits = c(0, 1), breaks = c(0, 0.5, 1), position = "top",
                     expand = expansion(mult = c(0, 0.02))) +
  geom_hline(yintercept = breaks_y$y, color = "#4d565f", linewidth = 0.4) +
  labs(x = "Communality", y = NULL) +
  theme_minimal(base_size = 12) +
  theme(axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(color = "#4d565f", size = 10))

p_heat + p_comm + plot_layout(widths = c(4, 1.3))

For a smaller set of variables, psych::fa.diagram(fa_fit) draws a path diagram from each factor to its indicators. It shows only each variable’s largest loading, though, so it hides exactly the cross-loadings the heatmap reveals.

A clean loadings table

Join the loadings back to the data dictionary, add the communality and uniqueness, and blank out loadings below 0.30 so the structure stands out.

L <- unclass(fa_fit$loadings)[, names(factor_names)]
colnames(L) <- factor_names

loadings_out <- as.data.frame(round(L, 2)) %>%
  rownames_to_column("LABEL") %>%
  mutate(
    communality = round(fa_fit$communality[LABEL], 2),
    uniqueness  = round(fa_fit$uniquenesses[LABEL], 2)
  ) %>%
  left_join(dd, by = "LABEL") %>%
  relocate(LABEL, VARIABLE)

readr::write_csv(loadings_out, file.path(out_dir, "factor_loadings.csv"))

loadings_out %>%
  mutate(across(all_of(unname(factor_names)),
                ~ ifelse(abs(.x) < 0.30, "", sprintf("%.2f", .x)))) %>%
  arrange(desc(communality)) %>%
  pander(justify = "left", split.table = Inf)

LABEL	VARIABLE	ses	stability	immigrant	aging	communality	uniqueness
polang12	Percent speaking other language at home, age 5 plus			0.98		1	0
p60up12	Percent 60 and older, total		0.33		0.92	1	0
pcol12	Percent with 4-year college degree or more	0.94				0.91	0.09
pmar12	Percent currently married, not separated	0.32	0.88			0.88	0.12
pfb12	Percent foreign born			0.86		0.85	0.15
hinc12	Median HH income, total	0.75	0.45			0.83	0.17
incpc12	Per capita income	0.86				0.81	0.19
phs12	Percent with high school degree or less	-0.82				0.77	0.23
p75up12	Percent 75 and older, total				0.85	0.76	0.24
pprof12	Percent professional employees	0.84				0.75	0.25
pown12	Percent owner-occupied units		0.84			0.75	0.25
phisp12	Percent Hispanic			0.84		0.74	0.26
pnhwht12	Percent white, non-Hispanic		0.56	-0.53		0.66	0.34
pmulti12	Percent multi-family units		-0.74			0.66	0.34
mhmval12	Median home value	0.76				0.65	0.35
ppov12	Percent in poverty, total	-0.51	-0.55			0.57	0.43
mrent12	Median rent	0.66		0.31		0.55	0.45
pfhh12	Percent female-headed families with children	-0.41	-0.55			0.49	0.51
pwds12	Percent widowed, divorced and separated	-0.35			0.53	0.44	0.56
pvet12	Percent veteran		0.37	-0.35	0.32	0.37	0.63
pnhblk12	Percent black, non-Hispanic		-0.52			0.36	0.64
psemp12	Percent self-employed	0.40	0.37			0.33	0.67
p18und12	Percent 17 and under, total				-0.41	0.33	0.67
punemp12	Percent unemployed	-0.40	-0.33			0.29	0.71
pflabf12	Percent female labor force participation				-0.42	0.26	0.74
p10yrs12	Percent HH in neighborhood 10 years or less		-0.31			0.17	0.83
pvac12	Percent vacant units				0.35	0.16	0.84
pmanuf12	Percent manufacturing employees					0.15	0.85
p30old12	Percent structures more than 30 years old					0.12	0.88
pntv12	Percent Native American race					0.01	0.99

The low end of the communality column is informative too. Percent Native American, percent of structures more than 30 years old, and percent manufacturing employees are barely explained by any of the four factors. They measure something else.

Top-loading variables for each factor

top_k <- 6

top_by_factor <- loadings_out %>%
  select(LABEL, VARIABLE, all_of(unname(factor_names))) %>%
  pivot_longer(all_of(unname(factor_names)), names_to = "factor", values_to = "loading") %>%
  group_by(factor) %>%
  slice_max(abs(loading), n = top_k, with_ties = FALSE) %>%
  ungroup() %>%
  mutate(factor = factor(factor, levels = unname(factor_names))) %>%
  arrange(factor, desc(abs(loading)))

readr::write_csv(top_by_factor, file.path(out_dir, "top_loadings_by_factor.csv"))

for (f in levels(top_by_factor$factor)) {
  rows <- filter(top_by_factor, factor == f)
  cat("\n#### ", f, "\n\n", sep = "")
  cat(paste0("- ", rows$VARIABLE, " (", sprintf("%.2f", rows$loading), ")", collapse = "\n"), "\n")
}

ses

Percent with 4-year college degree or more (0.94)
Per capita income (0.86)
Percent professional employees (0.84)
Percent with high school degree or less (-0.82)
Median home value (0.76)
Median HH income, total (0.75)

stability

Percent currently married, not separated (0.88)
Percent owner-occupied units (0.84)
Percent multi-family units (-0.74)
Percent white, non-Hispanic (0.56)
Percent in poverty, total (-0.55)
Percent female-headed families with children (-0.55)

immigrant

Percent speaking other language at home, age 5 plus (0.98)
Percent foreign born (0.86)
Percent Hispanic (0.84)
Percent white, non-Hispanic (-0.53)
Percent veteran (-0.35)
Median rent (0.31)

aging

Percent 60 and older, total (0.92)
Percent 75 and older, total (0.85)
Percent widowed, divorced and separated (0.53)
Percent female labor force participation (-0.42)
Percent 17 and under, total (-0.41)
Percent vacant units (0.35)

8. Does the solution depend on the estimator?

Because of the Heywood warning, refit the model with maximum likelihood (fm = "ml") and compare. factor.congruence() measures how similar two sets of loadings are. Values above 0.95 mean the factors are essentially the same.

fa_ml <- psych::fa(X_z, nfactors = n_factors, fm = "ml", rotate = "varimax")

round(psych::factor.congruence(fa_fit, fa_ml), 2)

##       ML3   ML4   ML1   ML2
## PA1  1.00  0.27 -0.02 -0.14
## PA2  0.46  0.99 -0.37  0.08
## PA4 -0.04 -0.33  1.00 -0.26
## PA3 -0.15  0.05 -0.25  1.00

Each principal axis factor has a maximum likelihood twin with a congruence of about 0.99 or higher, so the four constructs are not an artifact of the estimator.

9. Extract the new features

Factor scores are the payoff: one score per factor for every neighborhood. They replace 30 overlapping variables with four composite measures.

scores <- as.data.frame(fa_fit$scores)[, names(factor_names)]
names(scores) <- paste0("factor_", factor_names)

head(round(scores, 2))

factor_ses	factor_stability	factor_immigrant	factor_aging
1.43	0.75	-0.54	0.01
0.64	0.48	-0.53	0.74
1.80	-2.30	0.26	0.79
0.11	-0.52	-0.62	0.21
-0.63	0.47	-0.81	0.59
-0.27	2.21	-0.51	-1.59

readr::write_csv(scores, file.path(out_dir, "factor_scores.csv"))

dat_with_scores <- bind_cols(dat, scores)
readr::write_csv(dat_with_scores, file.path(out_dir, "data_with_factor_scores.csv"))

Under varimax the underlying factors are uncorrelated. The estimated scores are close to uncorrelated but not exactly, because each score is an estimate that contains some error.

scores_cor <- round(cor(scores), 3)
scores_cor

##                  factor_ses factor_stability factor_immigrant factor_aging
## factor_ses            1.000           -0.002            0.000        0.013
## factor_stability     -0.002            1.000            0.026       -0.031
## factor_immigrant      0.000            0.026            1.000        0.007
## factor_aging          0.013           -0.031            0.007        1.000

readr::write_csv(
  as.data.frame(scores_cor) %>% rownames_to_column("factor"),
  file.path(out_dir, "factor_score_correlations.csv")
)

10. Weighted scores versus a simple index

How much does factoring buy you over the simple approach, averaging the z-scores of the items that define a construct? Take the socioeconomic status factor. Build a simple index from its six strongest variables, flipping the sign of the one that loads negatively, and compare it to the factor score.

ses_items <- top_by_factor %>% filter(factor == "ses")
ses_items %>% select(LABEL, VARIABLE, loading)

LABEL	VARIABLE	loading
pcol12	Percent with 4-year college degree or more	0.94
incpc12	Per capita income	0.86
pprof12	Percent professional employees	0.84
phs12	Percent with high school degree or less	-0.82
mhmval12	Median home value	0.76
hinc12	Median HH income, total	0.75

# Reverse-code negatively loading items so every item points the same direction
ses_z <- X_z[ses_items$LABEL]
ses_z <- sweep(ses_z, 2, sign(ses_items$loading), `*`)

ses_simple <- rowMeans(ses_z)

# Cronbach's alpha for the unit-weighted scale
ses_alpha <- psych::alpha(ses_z)
round(ses_alpha$total$raw_alpha, 2)

## [1] 0.94

# How similar are the simple index and the weighted factor score?
round(cor(ses_simple, scores$factor_ses), 3)

## [1] 0.957

plot(ses_simple, scores$factor_ses,
     pch = 19, col = gray(0.5, 0.4), bty = "n",
     xlab = "Simple index (mean of six z-scores)",
     ylab = "Factor score (weighted)",
     main = "Socioeconomic status: two ways to build the index")
abline(lm(scores$factor_ses ~ ses_simple), col = "#bd5416", lwd = 2)

The six items are highly reliable as a scale (α = 0.94), and the simple index and the factor score are correlated 0.96. When items are strong and loadings are similar, a simple average and a weighted composite tell nearly the same story. Weighting matters most when some items are much noisier than others.

Outputs

Knitting this file writes five files to factor_analysis_outputs/.

To interpret and label each construct:

factor_loadings.csv: all loadings with communalities, uniquenesses, and descriptions
top_loadings_by_factor.csv: the six strongest variables for each factor

The new composite measures:

factor_scores.csv: four factor scores for each neighborhood
data_with_factor_scores.csv: the original data with the scores appended

A check on independence:

factor_score_correlations.csv

Notes on refinement

This is a teaching example that uses a four-factor solution because the loadings are easy to interpret. For a real instrument:

Drop weak indicators. Remove variables with low item-level KMO or low communality, such as pntv12, then refit.
Resolve redundancy. Pairs like percent foreign born and percent speaking another language (r = 0.92), or percent 60 and older and percent 75 and older (r = 0.87), carry nearly the same information. Keep one of each pair, or combine them.
Drop strong cross-loaders if you want a cleaner simple structure. Median household income and percent currently married both load on two factors.
Consider an oblique rotation. With rotate = "oblimin" the factors are allowed to correlate, and the factor correlation matrix ($Phi) tells you whether independence was a reasonable assumption.
Use maximum likelihood (fm = "ml") when you want formal fit statistics. It assumes approximate multivariate normality.
Confirm the structure on new data with a confirmatory factor analysis before treating the factors as validated measures.