---
title: "Visualizing Factor Analysis"
author: "Michael Hallquist, PSYC 859"
date: last-modified
week: 12
download: "materials/downloads/w12_models/visualizing_factor_analysis.qmd"
format:
  html:
    code-fold: show
execute:
  warning: false
  message: false
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.width = 10, fig.height = 7, echo = TRUE)
```

# Overview

Factor analysis sits in a different corner of the statistical modeling world than regression, but the core logic of visualization is the same: we have a **model** that makes predictions about data, and we want to see how well those predictions match reality.

In regression, the model predicts values of an outcome variable. In factor analysis, the model predicts the **pattern of correlations** among a set of observed variables. The "predictions" are the correlations implied by the factor solution, and "misfit" shows up as residual correlations — observed relationships that the factors cannot explain.

This walkthrough covers:

1. **Pre-analysis visualization** — examining the correlation structure before fitting a model.
2. **Determining the number of factors** — scree plots and parallel analysis.
3. **Interpreting the factor solution** — loading heatmaps and path diagrams.
4. **Visualizing fit and misfit** — residual correlation matrices that reveal what the model misses.
5. **Using factor scores** — connecting latent variables back to observable data.

## Learning goals

After working through this document, you should be able to:

1. Create and interpret a correlation heatmap for a set of questionnaire items.
2. Use scree plots and parallel analysis to decide on the number of factors.
3. Build a loading heatmap in `ggplot2` that reveals simple structure.
4. Compute and visualize residual correlations to assess model fit.
5. Extract factor scores and use them in downstream analyses.

## Packages

```{r}
# Install any missing packages (uncomment if needed):
# install.packages(c(
#   "ggplot2", "dplyr", "tidyr", "psych",
#   "patchwork", "scales"
# ))

library(ggplot2)
library(dplyr)
library(tidyr)
library(tibble)
library(psych)
library(patchwork)
library(scales)
```

# Data: The Big Five Inventory

We use the `bfi` dataset from the `psych` package. This dataset contains responses from 2,800 individuals to 25 personality items, with 5 items per Big Five trait: Agreeableness (A1–A5), Conscientiousness (C1–C5), Extraversion (E1–E5), Neuroticism (N1–N5), and Openness (O1–O5).

The BFI is ideal for demonstrating factor analysis because the factor structure is well established and psychologically interpretable.

```{r}
data(bfi, package = "psychTools")

raw_items <- bfi %>%
  select(matches("^[ACENO][0-9]+$"))

trait_labels <- c(
  A = "Agreeableness",
  C = "Conscientiousness",
  E = "Extraversion",
  N = "Neuroticism",
  O = "Openness"
)

score_keys <- psychTools::bfi.keys
names(score_keys) <- unname(trait_labels[c("A", "C", "E", "N", "O")])

reverse_key <- unlist(score_keys)
reverse_multiplier <- ifelse(grepl("^-", reverse_key), -1, 1)
names(reverse_multiplier) <- sub("^-", "", reverse_key)

# 1. Reverse negative items upfront
# This aligns all items positively with their trait label.
bfi_items <- psych::reverse.code(
  keys = reverse_multiplier[names(raw_items)],
  items = raw_items,
  mini = 1,
  maxi = 6
) %>%
  as.data.frame() %>%
  setNames(names(raw_items)) %>%
  as_tibble()

# 2. Re-attach demographics
# Notice we DO NOT drop_na(). We preserve all 2,800 rows.
bfi_data <- bind_cols(
  bfi_items,
  bfi %>% select(gender, age, education) %>%
    mutate(gender_f = factor(gender, levels = 1:2, labels = c("Male", "Female")))
)

cat("Total observations:", nrow(bfi_data), "\n")
cat("Items:", ncol(bfi_items), "\n")
```

# Pre-analysis: Examining the correlation structure

Before we fit any factor model, we should look at the correlation matrix. If items within the same trait cluster together, there is reason to expect that a factor model can capture that structure.

## Correlation heatmap with `psych::corPlot()`

The `psych` package offers a quick way to visualize correlation matrices. Items are ordered by their position in the questionnaire, which corresponds to the Big Five groupings:

```{r}
#| fig-width: 10
#| fig-height: 10
corPlot(bfi_items, upper = FALSE, diag = FALSE, main = "Correlation matrix of 25 BFI items", 
        gr = colorRampPalette(c("blue", "white", "red")), scale = FALSE)

```

You should see blocks of stronger positive correlations (red) along the diagonal — these correspond to items within the same trait. Because we reverse-coded all negatively keyed items in the step above, there are no longer any dark blue (negative) blocks within traits. Off-diagonal correlations between traits should be weaker (though not zero, because the Big Five are not perfectly orthogonal).

## Custom ggplot heatmap

For more control over the display, we can build a correlation heatmap in `ggplot2`:

```{r}
#| fig-width: 10
#| fig-height: 10
cor_mat <- cor(bfi_items, use = "pairwise.complete.obs")

# Convert to long format for ggplot
cor_long <- cor_mat %>%
  as.data.frame() %>%
  mutate(item1 = rownames(.)) %>%
  pivot_longer(-item1, names_to = "item2", values_to = "r") %>%
  # Create ordered factors so items appear in questionnaire order

  mutate(
    item1 = factor(item1, levels = colnames(cor_mat)),
    item2 = factor(item2, levels = rev(colnames(cor_mat)))
  )

ggplot(cor_long, aes(x = item1, y = item2, fill = r)) +
  geom_tile(color = "white", linewidth = 0.3) +
  scale_fill_gradient2(
    low = "#2166ac", mid = "white", high = "#b2182b",
    midpoint = 0, limits = c(-1, 1),
    name = "r"
  ) +
  labs(
    title = "Correlation matrix of 25 BFI items",
    subtitle = "Blue = negative; Red = positive; Items grouped by trait",
    x = NULL, y = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank()
  ) +
  coord_fixed()
```

# Determining the number of factors

## Scree plot

A scree plot shows the eigenvalues of the correlation matrix in descending order. The idea is to look for an "elbow" — a point where eigenvalues drop off sharply. Factors before the elbow capture substantial shared variance; factors after it capture noise.

```{r}
eigen_vals <- eigen(cor_mat)$values

scree_df <- data.frame(
  factor = 1:length(eigen_vals),
  eigenvalue = eigen_vals
)

ggplot(scree_df, aes(x = factor, y = eigenvalue)) +
  geom_line(color = "steelblue", linewidth = 1) +
  geom_point(color = "steelblue", size = 3) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "gray50") +
  scale_x_continuous(breaks = 1:25) +
  labs(
    title = "Scree plot",
    subtitle = "Dashed line at eigenvalue = 1 (Kaiser criterion)",
    x = "Factor number",
    y = "Eigenvalue"
  ) +
  theme_minimal(base_size = 14)
```

The Kaiser criterion (eigenvalue > 1) and the scree plot often suggest slightly different numbers. Horn's parallel analysis provides a more rigorous approach.

## Parallel analysis

The standard scree plot has a major limitation: even a matrix of completely random noise will produce some factors with eigenvalues greater than 1 simply due to chance sampling fluctuations. 

**Parallel analysis** solves this by simulating completely random datasets that have the exact same number of participants and items as your real data. It calculates the eigenvalues for these random datasets, allowing us to ask: *Is the variance captured by our real factor larger than the variance we would expect just by chance?*

```{r}
#| fig-width: 10
#| fig-height: 7
fa_parallel <- fa.parallel(bfi_items, fa = "fa", fm = "minres", main = "Parallel analysis of BFI items")
```

### How to interpret the graph

The plot produced by `fa.parallel()` contains three types of lines:

1. **Actual Data (solid blue line, triangles):** These are the true eigenvalues extracted from your dataset. They represent the actual variance captured by each successive factor.
2. **FA Simulated Data (dashed red line):** These are the average eigenvalues extracted from many simulated datasets of pure random noise (where the true correlation between all items is exactly zero). 
3. **FA Resampled Data (dotted red line):** These are eigenvalues generated by randomly shuffling your actual data across participants. This preserves the individual distributions of your items while destroying any correlations between them.

**The decision rule:** You retain only the factors where the **solid blue line (Actual Data)** is *above* the red lines (Simulated and Resampled data). Once the blue line drops below the red lines, the remaining "factors" in your real data are capturing less variance than factors extracted from literal random noise.

In this plot, the blue line drops below the red lines after the 6th factor. However, parallel analysis can sometimes slightly over-extract.

## A consensus approach: MAP and VSS

No single test for the number of factors is perfect. As psychometricians often note, parallel analysis can slightly over-extract, while Velicer's **MAP (Minimum Average Partial)** test can sometimes under-extract. The **VSS (Very Simple Structure)** criterion assesses how well the data are fit by factors of varying complexity.

Rather than relying purely on parallel analysis, the most sensible approach is to look for consensus across multiple indices. The `psych::nfactors()` function runs all of these tests at once:

```{r}
# Run multiple criteria (MAP, VSS, and eigenvalue tests)
n_factors_check <- nfactors(bfi_items, n = 8, rotate = "oblimin", fm = "minres")
n_factors_check
```

### Interpreting the `nfactors` plot

When you run `nfactors()`, it generates a four-panel diagnostic plot summarizing the fit across different numbers of factors:

1. **Very Simple Structure (VSS):** Shows the VSS fit for solutions restricted to a complexity of 1 (each item loads on only one factor) or 2. You are looking for the peak—the number of factors where the VSS fit is *maximized*.
2. **Complexity:** Shows Hoffman's index of complexity. A lower value indicates a "cleaner", simpler structure where items don't cross-load heavily onto multiple factors.
3. **Empirical BIC:** The Bayesian Information Criterion balances model fit against complexity (penalizing for too many factors). You are looking for the *minimum* (lowest overall point) on the curve.
4. **Root Mean Residual (RMSR):** Shows the average magnitude of the residual correlations (the model deviations). Lower values mean the model reproduces the observed data more accurately.

By examining the output and the plots, you can see if MAP, empirical BIC, and VSS converge with parallel analysis and traditional scree rules. If they point to different numbers (e.g., MAP suggests 5, parallel suggests 6), you must use theoretical knowledge and model residuals to make the final choice. Given the strong theoretical support for the Big Five, and since the BFI dataset was designed explicitly around 5 traits, we will proceed with the 5-factor solution.

# Fitting the factor model

We use `psych::fa()` with oblimin rotation (which allows factors to correlate — appropriate for personality traits, which are not independent).

```{r}
fa_5 <- fa(bfi_items, nfactors = 5, rotate = "oblimin", fm = "minres")
print(fa_5, cut = 0.3, sort = TRUE)
```

## Path diagram

`psych::fa.diagram()` produces a path diagram showing which items load on which factors:

```{r}
#| fig-width: 10
#| fig-height: 8
fa.diagram(fa_5, main = "Factor structure: 5-factor oblimin solution")
```

## Loading heatmap in ggplot

A custom heatmap of factor loadings gives us more control and a cleaner display. We want items sorted by their primary loading so that simple structure is visually apparent.

```{r}
#| fig-width: 8
#| fig-height: 12

# Extract loadings matrix
loadings_mat <- as.data.frame(unclass(fa_5$loadings))
loadings_mat$item <- rownames(loadings_mat)

# Determine each item's primary factor (highest absolute loading)
loadings_mat$primary_factor <- apply(abs(loadings_mat %>% select(-item)), 1, which.max)

loadings_mat$max_loading <- apply(abs(loadings_mat %>% select(-item, -primary_factor)), 1, max)

# Sort items: first by primary factor, then by loading magnitude within factor
loadings_mat <- loadings_mat %>%
  arrange(primary_factor, desc(max_loading))

# Convert to long format
loadings_long <- loadings_mat %>%
  select(-primary_factor, -max_loading) %>%
  pivot_longer(-item, names_to = "factor", values_to = "loading") %>%
  mutate(item = factor(item, levels = rev(loadings_mat$item)))

# Rename factors for interpretability
factor_labels <- c(
  "MR1" = "F1", "MR2" = "F2", "MR3" = "F3",
  "MR4" = "F4", "MR5" = "F5"
)
loadings_long$factor <- factor_labels[loadings_long$factor]

ggplot(loadings_long, aes(x = factor, y = item, fill = loading)) +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(
    aes(label = ifelse(abs(loading) >= 0.3, sprintf("%.2f", loading), "")),
    size = 3, color = "black"
  ) +
  scale_fill_gradient2(
    low = "#2166ac", mid = "white", high = "#b2182b",
    midpoint = 0, limits = c(-1, 1),
    name = "Loading"
  ) +
  labs(
    title = "Factor loading heatmap",
    subtitle = "Items sorted by primary factor; values shown for |loading| ≥ .30",
    x = "Factor",
    y = "Item"
  ) +
  theme_minimal(base_size = 13) +
  theme(panel.grid = element_blank())
```

A well-fitting factor model will show **simple structure**: each item loads strongly on one factor and weakly on the others. Items with strong cross-loadings suggest measurement problems or factors that are not cleanly separated.

# Visualizing fit and misfit

## The logic of residual correlations

A factor model's fundamental prediction is about the correlation matrix. The model says:

$$\hat{\mathbf{R}} = \mathbf{\Lambda} \mathbf{\Phi} \mathbf{\Lambda}' + \mathbf{\Psi}$$

where $\hat{\mathbf{R}}$ is the model-implied correlation matrix, $\mathbf{\Lambda}$ is the loading matrix, $\mathbf{\Phi}$ is the factor correlation matrix, and $\mathbf{\Psi}$ is the uniqueness (residual variance) diagonal.

The **residual correlation matrix** is simply:

$$\mathbf{R}_{residual} = \mathbf{R}_{observed} - \hat{\mathbf{R}}$$

Large residual correlations flag item pairs whose relationship the factor model does not explain. These are the places where the model is *wrong*.

## Computing and plotting residual correlations

```{r}
#| fig-width: 10
#| fig-height: 10

# Extract residual correlation matrix from the fa object
resid_mat <- fa_5$residual

# Zero out the diagonal (residual variances are not interesting here)
diag(resid_mat) <- NA

# Convert to long format
resid_long <- resid_mat %>%
  as.data.frame() %>%
  mutate(item1 = rownames(.)) %>%
  pivot_longer(-item1, names_to = "item2", values_to = "residual") %>%
  mutate(
    item1 = factor(item1, levels = colnames(resid_mat)),
    item2 = factor(item2, levels = rev(colnames(resid_mat)))
  )

ggplot(resid_long, aes(x = item1, y = item2, fill = residual)) +
  geom_tile(color = "white", linewidth = 0.3) +
  scale_fill_gradient2(
    low = "#2166ac", mid = "white", high = "#b2182b",
    midpoint = 0, limits = c(-0.15, 0.15), oob = squish,
    name = "Residual\ncorrelation"
  ) +
  labs(
    title = "Residual correlations: What the 5-factor model misses",
    subtitle = "Large colored cells = item pairs not well explained by 5 factors",
    x = NULL, y = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank()
  ) +
  coord_fixed()
```

In a well-fitting model, the residual correlation matrix should be mostly pale (near zero). Pockets of strong colour indicate item pairs whose relationship exceeds what the factors predict.

## Fit indices

The numeric fit indices connect directly to the visual residuals:

```{r}
cat("Root Mean Square Residual (RMSR):", round(fa_5$rms, 4), "\n")
cat("RMSEA:", round(fa_5$RMSEA["RMSEA"], 4),
    " [90% CI:", round(fa_5$RMSEA["lower"], 4),
    ",", round(fa_5$RMSEA["upper"], 4), "]\n")
cat("TLI:", round(fa_5$TLI, 4), "\n")
cat("BIC:", round(fa_5$BIC, 2), "\n")
```

Conventional guidelines (which should be interpreted cautiously):

- **RMSR**: Should be close to the expected value for your sample size (printed by `fa()`). Lower is better.
- **RMSEA**: < .05 is "good," .05–.08 is "acceptable." Higher indicates misfit.
- **TLI**: > .90 is "acceptable," > .95 is "good."

# Comparing factor solutions

One powerful diagnostic is to compare the residual correlations across different numbers of factors. If adding a factor eliminates a block of residual correlations, it was capturing real covariance.

```{r}
fa_4 <- fa(bfi_items, nfactors = 4, rotate = "oblimin", fm = "minres")
fa_6 <- fa(bfi_items, nfactors = 6, rotate = "oblimin", fm = "minres")
```

## Side-by-side comparison of fit

```{r}
comparison <- data.frame(
  Factors = c(4, 5, 6),
  RMSR = c(fa_4$rms, fa_5$rms, fa_6$rms),
  RMSEA = c(fa_4$RMSEA["RMSEA"], fa_5$RMSEA["RMSEA"], fa_6$RMSEA["RMSEA"]),
  TLI = c(fa_4$TLI, fa_5$TLI, fa_6$TLI),
  BIC = c(fa_4$BIC, fa_5$BIC, fa_6$BIC)
)

comparison %>%
  mutate(across(RMSR:TLI, ~ round(.x, 4)), BIC = round(BIC, 1))
```

## Residual correlation comparison

```{r}
#| fig-width: 14
#| fig-height: 5

# Helper function to create a residual heatmap
make_resid_heatmap <- function(fa_obj, title_label) {
  resid_m <- fa_obj$residual
  diag(resid_m) <- NA

  resid_m %>%
    as.data.frame() %>%
    mutate(item1 = rownames(.)) %>%
    pivot_longer(-item1, names_to = "item2", values_to = "residual") %>%
    mutate(
      item1 = factor(item1, levels = colnames(resid_m)),
      item2 = factor(item2, levels = rev(colnames(resid_m)))
    ) %>%
    ggplot(aes(x = item1, y = item2, fill = residual)) +
    geom_tile(color = "white", linewidth = 0.2) +
    scale_fill_gradient2(
      low = "#2166ac", mid = "white", high = "#b2182b",
      midpoint = 0, limits = c(-0.15, 0.15), oob = squish
    ) +
    labs(title = title_label, x = NULL, y = NULL) +
    theme_minimal(base_size = 9) +
    theme(
      axis.text.x = element_text(angle = 90, hjust = 1, size = 7),
      axis.text.y = element_text(size = 7),
      panel.grid = element_blank(),
      legend.position = "none"
    ) +
    coord_fixed()
}

p4 <- make_resid_heatmap(fa_4, "4 factors")
p5 <- make_resid_heatmap(fa_5, "5 factors")
p6 <- make_resid_heatmap(fa_6, "6 factors")

p4 + p5 + p6 +
  plot_annotation(
    title = "Residual correlations by number of factors",
    subtitle = "More factors → smaller residuals (but diminishing returns)"
  )
```

As the number of factors increases, the residual correlation matrix should become increasingly pale. But there are diminishing returns: if adding a 6th factor only eliminates a tiny pocket of residual color, it may not be worth the added complexity.

# Factor scores

Once we are satisfied with the factor solution, we can extract **factor scores** — each person's estimated standing on each latent factor. Alternatively, we can use classical test theory to compute simple sum or mean scores across the items that make up each factor.

## 1. Classical unit-weighted scores (`scoreItems`)

A common alternative to EFA-derived factor scores is to simply average the raw items that belong to each trait. The `psych::scoreItems()` function does this automatically. Crucially, if we specify `impute = "median"`, it will handle missing data by substituting the median response for that person's missing item, allowing us to keep all 2,800 participants.

```{r}
# Since we already reverse-scored the items upfront, our keys just point 
# to the existing positive columns for each trait.
scale_keys <- list(
  Agreeableness = paste0("A", 1:5),
  Conscientiousness = paste0("C", 1:5),
  Extraversion = paste0("E", 1:5),
  Neuroticism = paste0("N", 1:5),
  Openness = paste0("O", 1:5)
)

classic_scores <- scoreItems(
  keys = scale_keys, 
  items = bfi_items, 
  impute = "median"
)

# Extract the scores matrix into a data frame
scores_classic <- as.data.frame(classic_scores$scores) %>% 
  rename_with(~ paste0(.x, "_classic"))
```

## 2. EFA factor scores (`factor.scores`)

We can also extract the mathematically derived factor scores from the EFA model. We use the `tenBerge` method because it guarantees that the extracted scores will have the exact same correlations as the latent factors. We also pass `impute = "median"` so that individuals with missing item data are retained rather than dropped via listwise deletion.

```{r}
fs_results <- factor.scores(bfi_items, fa_5, method = "tenBerge", impute = "median")
scores_efa <- as.data.frame(fs_results$scores)

# Rename the columns to match their psychological traits based on item loadings
# (Check the loading matrix: MR2=Neuroticism, MR1=Extraversion, etc.)
colnames(scores_efa) <- c("Neuroticism_efa", "Extraversion_efa", "Conscientiousness_efa", "Agreeableness_efa", "Openness_efa")

# Bind both score types back to our original dataset
bfi_scored <- bind_cols(bfi_data, scores_classic, scores_efa)
```

> **Identifying the factors:** Notice that `fa()` output calls the factors MR1, MR2, etc. These numbers are *arbitrary* and depend on the extraction math, not psychological theory. To figure out "Which factor is which?", look at the heatmaps or the path diagram to see which items load highest. For example, if the Agreeableness items (A1-A5) all load strongly on MR2, then MR2 *is* the Agreeableness factor. It's a best practice to rename factors so they are interpretable.

## Factor score distributions by gender

```{r}
#| fig-width: 10
#| fig-height: 8
scores_long <- bfi_scored %>%
  filter(!is.na(gender_f)) %>%
  pivot_longer(ends_with("_efa"), names_to = "factor", values_to = "score") %>%
  mutate(factor = sub("_efa", "", factor))

ggplot(scores_long, aes(x = score, fill = gender_f)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~factor, scales = "free_y", ncol = 1) +
  scale_fill_manual(
    values = c("Male" = "#2c7fb8", "Female" = "#d95f02"),
    name = "Gender"
  ) +
  labs(
    title = "Factor score distributions by gender",
    subtitle = "Each panel is one Big Five factor",
    x = "Factor score",
    y = "Density"
  ) +
  theme_minimal(base_size = 13) +
  theme(strip.text = element_text(face = "bold"))
```

## Factor scores as regression predictors

Factor scores can be used as predictors or outcomes in subsequent analyses, connecting the latent measurement model back to the regression framework covered in the first walkthrough.

```{r}
# Example: does neuroticism predict agreeableness?
m_fs <- lm(Agreeableness_efa ~ Neuroticism_efa + age + gender_f, data = bfi_scored)
summary(m_fs)
```

```{r}
# Quick visualization
ggplot(bfi_scored %>% filter(!is.na(gender_f)), aes(x = Neuroticism_efa, y = Agreeableness_efa, color = gender_f)) +
  geom_point(alpha = 0.15, size = 1) +
  geom_smooth(method = "lm", se = TRUE) +
  scale_color_manual(values = c("Male" = "#2c7fb8", "Female" = "#d95f02")) +
  labs(
    title = "Factor scores: Neuroticism predicting Agreeableness",
    subtitle = "Controlling for gender (separate regression lines)",
    x = "Neuroticism score",
    y = "Agreeableness score",
    color = "Gender"
  ) +
  theme_minimal(base_size = 14)
```

# Take-home summary

## The factor analysis visualization workflow

1. **Before fitting**: Examine the correlation matrix with a heatmap. Look for blocks of correlated items that suggest latent factors. Use parallel analysis to determine the number of factors.
2. **After fitting**: Visualize the loading matrix as a heatmap. Look for simple structure: each item loading strongly on one factor. Use `fa.diagram()` for a quick structural overview.
3. **Diagnosing misfit**: Compute the residual correlation matrix (observed minus model-implied). Large residual correlations reveal item pairs that the model fails to explain. Compare residual matrices across different numbers of factors.
4. **Using the results**: Extract factor scores and visualize their distributions. Factor scores can then be used in regression models, group comparisons, or other analyses.

## Tool reference

| Task | Function | Package |
|---|---|---|
| Correlation heatmap | `corPlot()` or custom `ggplot` with `geom_tile()` | psych / ggplot2 |
| Scree plot + parallel analysis | `fa.parallel()` | psych |
| Exploratory factor analysis | `fa()` | psych |
| Path diagram | `fa.diagram()` | psych |
| Loading heatmap | Custom `ggplot` with `geom_tile()` + `geom_text()` | ggplot2 |
| Residual correlation heatmap | Custom `ggplot` using `fa_obj$residual` | ggplot2 |
| Factor score extraction | `fa_obj$scores` | psych |
| Factor score distributions | `geom_density()` / `geom_violin()` | ggplot2 |

## Further reading

- Revelle, W. (2024). *psych: Procedures for personality and psychological research*. <https://personality-project.org/r/psych/>
- Revelle, W. (2024). *An introduction to psychometric theory with applications in R*. <http://personality-project.org/r/book/>
- The "Big Five" personality structure was originally described by Goldberg (1990) and the BFI items were developed by John, Donahue, & Kentle (1991).
