---
title: "Chapter 17: Exploratory Factor Analysis"
author: "Andy Field"
date: "2025-06-12"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
    code-overflow: wrap
execute: 
    message: false
    warning: false
---

## Setup

```{r}
#| label: setup

library(easystats)
library(tidyverse)


# load data

raq_tib <- here::here("data/raq.csv") |> 
  read_csv()

# or load from the discovr package
# raq_tib <- discovr::raq

```

## Prepare data

First we strip out items we don't need and reverse score negatively phrased items.

```{r}
#| label: raq_prep

# retain only items
raq_items <- raq_tib |> 
  select(-id)

# reverse score item 3
raq_items <- raq_items |>
  reverse(select = "raq_03", range = c(1, 5))
```

With more experience you'd combine these steps into a single pipe:

```{r}
#| label: data_prep_pipe
#| eval: false

raq_items <- here::here("data/raq.csv") |> 
  read_csv() |> 
  select(-id) |> 
  reverse(select = "raq_03", range = c(1, 5))

```


## Get a feel for the data

```{r}
#| label: cors

# compute polychoric correlations
raq_poly <- correlation(raq_items, method = "polychoric")
raq_poly_mtx <- as.matrix(raq_poly) # store correlation matrix

# display them in a matrix layout
raq_poly |>
  summary() |> 
  display(stars = FALSE)

# plot them

raq_poly |> 
  summary() |> 
  plot(stars = FALSE) +
  theme_minimal()
```

## Check data

```{r}
#| label: kmo

# get the KMO and Bartlett stats
check_factorstructure(raq_items)

# Determinant
raq_poly |> 
  as.matrix() |> 
  det()
```

## Distibrutions of items

```{r}
#| label: item_hists

raq_tidy_tib <- raq_items |>
  pivot_longer(
    cols = raq_01:raq_23,
    names_to = "Item",
    values_to = "Response"
  ) |> 
  mutate(
    Item = gsub("raq_", "RAQ ", Item)
  )

ggplot(raq_tidy_tib, aes(Response)) +
  geom_histogram(binwidth = 1, fill = "#CC6677", colour = "#882255", alpha = 0.5) +
  labs(y = "Frequency") +
  facet_wrap(~ Item, ncol = 6) +
  theme_minimal()


```

## Number of factors to extract

### Parallel analysis

```{r}
#| label: parallel

# from the data (slow)
# psych::fa.parallel(x = raq_items, cor = "poly")

# from the correlation matrix
psych::fa.parallel(x = raq_poly_mtx, n.obs = 2571, fa = "fa")
psych::fa.parallel(x = raq_poly_mtx, n.obs = 2571, fa = "pc")
```


### Method of agreement

```{r}
#| label: n_factors

raq_nfac <- n_factors(x = raq_items, correlation_matrix = raq_poly_mtx)

raq_nfac # written summary of results
display(raq_nfac) # table of methods/results
plot(raq_nfac) + theme_minimal() # plot of results
```

## Factor analysis
### Fit the model

```{r}
#| label: efa

# from raw scores:
raq_fa <- factor_analysis(x = raq_items,
                    n = 4,
                    scores = "tenBerge",
                    cor = "poly")
```

Alternatively, use the correlation matrix from earlier:

```{r}
#| label: efa_mat
#| eval: false

raq_fa <- factor_analysis(x = raq_poly_mtx,
                          n.obs = 2571,
                          n = 4)
```


### Check the model

```{r}
#| label: model_fit


model_performance(raq_fa) |> 
  display()

check_residuals(raq_fa)
check_normality(raq_fa) |> plot(detrend = FALSE)
```

:::{.callout-note}
## Bonus content!

The following code is for pedagogic purposes only.
:::

The following code explores what residuals are in factor analysis.

```{r}
#| label: explain_resids


# View a portion of the original r matrix
poly_mtx_mini <- raq_poly_mtx[1:5, 1:5] 
poly_mtx_mini |> 
  round(3)

# View a portion of the reproduced r matrix
repro_mtx <- get_model(raq_fa, element = "model") # extract reproduced correlation matrix
repro_mtx_mini <- repro_mtx[1:5, 1:5] # select only first 5 rows/columns
repro_mtx_mini |> 
  round(3)

# View the difference:
raq_resids <- poly_mtx_mini-repro_mtx_mini
raq_resids |> 
  round(3)

# get the residuals
fa_resids <- get_residuals(raq_fa)

# explaining what the RMSR is
fa_resids^2 |> 
  mean() |> 
  sqrt()
```

### Model summary

```{r}
#| label: efa_summary

summary(raq_fa) |> 
  display(digits = 2)
```



### Interpret the factors

```{r}
#| label: efa_load

raq_fa |> 
  display()
```

Suppress loadings below


```{r}
#| label: efa_reload

raq_fa |> 
  display(threshold = 0.2, sort = TRUE)
```

### Factor scores

```{r}
#| label: scores

factor_scores(raq_fa) |> 
  display()
```


## Omega

```{r}
#| label: omg
 
raq_omg <- item_omega(raq_items,
                      n = 4,
                      poly_cor = TRUE)
```


### Factor loadings

```{r}
#| label: omg_load

model_parameters(raq_omg) |> 
  display(threshold = 0.2)
```

### Reliability coefficients

```{r}
#| label: omg_sum

summary(raq_omg) |> 
  display()
```


### Fit statistics

```{r}
#| label: omg_fit

model_performance(raq_omg) |> 
  display()
```

## Cronbach's alpha

```{r}
#| label: alpha

item_to_factor <- c(3, 2, 3, 3, 3, 1, 1, 4, 2, 1, 4, 3, 1, 1, 1, 3, 4, 1, 2, 3, 3, 2, 2)

check_itemscale(x = raq_items, factor_index = item_to_factor) |> 
  display()

```
