Factor Analysis

Author

Riccardo Volpato

# OPTION 1: Simulate data
source("simulate_data.R")

efa_data <- simulate_data(
  sample_path = "data/TPSai+-+Pilot+Tester_27+April+2026_13.02.csv", 
  factors_path =  "data/factors_v1.1.csv",
  seed = seed,
  n = 250)
cols_tpsai <- names(efa_data)[c(2:40)]

Exploratory Factor Analysis

Brief Explainer

Exploratory factor analysis (EFA) is a statistical technique for identifying hidden structure in a set of items. The core idea is simple: if many items correlate with each other, they are probably all measuring some common underlying thing, a latent factor. EFA works by finding the smallest number of latent factors that can reproduce the observed correlations among items. It does this by asking: if these factors exist, what pattern of correlations between items would they produce? The solution that best matches the actual observed correlations, leaving the smallest residuals, is retained. Because EFA makes no assumptions about which items belong to which factor, it is well suited to the early stages of scale development, where the goal is to let the data reveal the underlying structure rather than confirm a predetermined one.

We use EFA to assess whether our theoretical factor structure is confirmed by the empirical factor structure in the data.

Factorability

Before applying EFA, we examine the Kaiser-Meyer-Olkin (KMO) index which asks whether the correlations among items are sufficiently compact and patterned to yield reliable factors. Values above .80 indicate good sampling adequacy and values above .90 are excellent. We treat KMO as our primary indicator of factorability.

Thorughout our EFA analysis, we use polychoric correlation matrix rather than the standard Pearson correlation matrix, as polychoric correlations are more appropriate for ordinal Likert data (Boateng et al. 2018).

R <- polychoric(efa_data %>% select(all_of(cols_tpsai)))$rho

KMO(R)$MSA

[1] 0.89585

Parallel Analysis

Brief Explainer

A key decision in EFA is how many factors to extract. Extracting too few factors leaves meaningful structure unexplained; extracting too many produces factors that are essentially fitting noise in the data. Parallel analysis addresses this by comparing each factor’s eigenvalue from the real correlation matrix to the corresponding eigenvalue from random correlation matrices of the same dimensions. Each factor is associated with an eigenvalue. Mathematically, an eigenvalue summarises how much of the shared variance among items a factor captures. When EFA extracts a factor, it finds a weighted combination of items that explains as much of the common variance as possible, the eigenvalue is the sum of the squared factor loadings for that factor across all items. A large eigenvalue therefore means many items have high loadings on that factor. Eigenvalues are computed using MINRES estimation, which minimises the residual correlations between the observed and reproduced correlation matrices and makes no distributional assumptions. In parallel analysis, factors are ordered by their eigenvalues from largest to smallest, with the first eigenvalue belonging to the first factor, the one explaining the most variance of the real correlation matrix, the second to the second factor, and so on. For each factor in turn, if the eigenvalue from the real data exceeds the average eigenvalue from the random matrices, meaning the factor is accounting for more variance than a factor extracted from pure noise would, that factor is considered to reflect genuine underlying structure in the data and is worth retaining. The moment the real eigenvalue drops below the random one, which is also what the scree plot visualises, with the blue line (real data) dropping below the red line (random data), that factor and all subsequent ones are considered noise and discarded.

fa <- fa.parallel(R, n.obs = nrow(efa_data), fm = "minres", fa = "fa")

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

Factor Estimation

Brief Explainer

Once the number of factors is determined, we extract them using MINRES estimation with oblimin rotation.

MINRES (minimum residual) estimation finds the factor solution that minimises the residuals, the differences between the observed correlation matrix and the one reproduced by the factor model. It makes no distributional assumptions, making it well suited to ordinal Likert data compared to maximum likelihood estimation, which assumes multivariate normality. It is the same method we used to compute eigenvalues and determine the number of factors.

Rotation is applied after initial factor extraction to make the solution more interpretable. Without rotation, factors are mathematically defined but often difficult to interpret because most items load moderately onto most factors. Rotation redistributes the variance among factors to achieve a simpler, cleaner pattern where each item loads strongly onto one factor and weakly onto others.

We use oblimin, which is an oblique rotation method meaning factors are allowed to correlate with one another. This is theoretically appropriate here because items within each factor are expected to measure related but distinct aspects of the same broader construct. An orthogonal rotation such as varimax would force factors to be uncorrelated, which is a stronger and less realistic assumption in psychological measurement. Importantly, if the factors turn out to be empirically uncorrelated, oblimin will recover the same solution as an orthogonal rotation, so nothing is lost by choosing oblique rotation.

efa <- fa(
  r=R, 
  nfactors = fa$nfact, # can ither user `fa$nfact` or visual 'elbow' criteria from plot
  n.obs = nrow(efa_data),
  fm = "minres",
  rotate = "oblimin"
)

Before exploring item loadings, we can visualise which items are assigned to which factor. This is also what we can use to export an EFA model for CFA.

efa_to_cfa(efa)

# Latent variables
MR1 =~ COMP_02 + COMP_03 + COMP_05 + CORP_03 + EMOT_04 + PERS_04 + PERS_05 + RISK_03 + SAFE_02
MR3 =~ COMP_01 + COMP_04 + CRIT_03 + EMOT_01 + EMOT_03 + PERS_01 + PERS_02 + PERS_03 + RISK_02 + RISK_04 + SAFE_03 + SAFE_05
MR2 =~ CORP_01 + CORP_04 + CORP_05 + HUMN_04 + HUMN_05 + RISK_05
MR7 =~ EMOT_02 + EMOT_05 + SAFE_01 + SAFE_04
MR5 =~ CRIT_01 + CRIT_04 + HUMN_01
MR6 =~ CRIT_05 + HUMN_02 + HUMN_03 + RISK_01
MR4 =~ CRIT_02

Factor Loadings

plot_loadings <- function(efa, order = c("theoretical", "empirical")) {
  
  order <- match.arg(order)
  
  loads <- efa$loadings |>
    unclass() |>
    as.data.frame() |>
    rownames_to_column("item")
  
  factor_cols <- setdiff(names(loads), "item")
  
  loads <- loads |>
    rowwise() |>
    mutate(
      primary_factor = factor_cols[which.max(abs(c_across(all_of(factor_cols))))],
      max_loading    = max(abs(c_across(all_of(factor_cols))))
    ) |>
    ungroup()
  
  ordered_items <- if (order == "empirical") {
    loads |>
      arrange(primary_factor, desc(max_loading)) |>
      pull(item)
  } else {
    loads |>
      pull(item)
  }
  
  loads |>
    pivot_longer(cols = -c(item, primary_factor, max_loading), names_to = "factor", values_to = "loading") |>
    mutate(
      item = factor(item, levels = ordered_items),
      label = ifelse(abs(loading) >= 0.3, as.character(round(loading, 2)), "-")
    ) |>
    ggplot(aes(x = factor, y = item, fill = loading)) +
    geom_tile() +
    geom_text(aes(label  = label, colour = abs(loading) > 0.5)) +
    scale_colour_manual(values = c("TRUE" = "white", "FALSE" = "black"), guide = "none") +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
    scale_y_discrete(limits = rev) +
    theme_minimal() +
    labs(x = "", y = "", fill = "")
}

We visualise how the items load across factors through a heat-map. First, are interested to see whether the empirical factor structure suggested by EFA matches our theoretical factor structure. Thus we plot the EFA results ordered by theoretical factors. In our plot, we hide the value of loadings below 0.3 to avoid visual cluttering.

plot_loadings(efa, order = "theoretical")

Since the empirical factor structure of EFA seems largely unrelated to our theoretical factor structure, we then visualise the factor structure with items ordered by empirical loadings.

plot_loadings(efa, order = "empirical")

An ideal EFA factors solution is robust and simple, that is when loadings have high values in one factor only (Buchanan, Valentine, and Schulenberg 2014). More specifically, in an good/interpretable solution each factor has at least 3 items and there are no “bad” items such as weak loaders that have a low primary loading (< 0.5) or cross-loaders that have small difference between primary and secondary loading (< 0.15).

flag_items <- function(efa, min_primary = 0.5, min_diff = 0.15) {
  
  efa$loadings |>
    unclass() |>
    as.data.frame() |>
    rownames_to_column("item") |>
    rowwise() |>
    mutate(
      abs_loads      = list(sort(abs(c_across(-item)), decreasing = TRUE)),
      primary        = abs_loads[[1]],
      secondary      = abs_loads[[2]],
      diff           = primary - secondary,
      low_primary    = primary < min_primary,
      low_diff       = diff < min_diff,
      flagged        = low_primary | low_diff
    ) |>
    ungroup() |>
    filter(flagged) |>
    select(item, primary, secondary, diff, low_primary, low_diff) |>
    knitr::kable(digits = 2)
}

flag_items(efa)

item	primary	secondary	diff	low_primary	low_diff
COMP_01	0.36	0.35	0.00	TRUE	TRUE
COMP_03	0.49	0.19	0.30	TRUE	FALSE
COMP_04	0.40	0.27	0.13	TRUE	TRUE
CORP_01	0.45	0.19	0.26	TRUE	FALSE
CORP_03	0.32	0.21	0.11	TRUE	TRUE
CRIT_01	0.46	0.29	0.18	TRUE	FALSE
CRIT_03	0.25	0.21	0.04	TRUE	TRUE
CRIT_04	0.28	0.27	0.01	TRUE	TRUE
EMOT_01	0.37	0.28	0.08	TRUE	TRUE
EMOT_03	0.45	0.23	0.21	TRUE	FALSE
EMOT_04	0.43	0.38	0.06	TRUE	TRUE
EMOT_05	0.47	0.40	0.07	TRUE	TRUE
HUMN_02	0.41	0.24	0.18	TRUE	FALSE
HUMN_03	0.45	0.32	0.13	TRUE	TRUE
HUMN_05	0.43	0.35	0.08	TRUE	TRUE
RISK_01	0.30	0.29	0.01	TRUE	TRUE
RISK_02	0.42	0.38	0.04	TRUE	TRUE
RISK_05	0.50	0.39	0.11	FALSE	TRUE
SAFE_02	0.47	0.31	0.16	TRUE	FALSE

Factor Diagnostics

In order to understand why our theoretical factor structure does not match the empirical factor structure suggested by EFA, we can inspect the structure of correlations in our data, including factor-factor, item-factor, and cross-factor inter-item correlations. By examining these we can search for a different factor structure that might have more empirical support from the collected data and revise our factor structure and its items accordingly.

Factor correlations

We report Spearman correlations between the theoretical factors to examine their inter-relationships. Spearman correlations are used as a rank-based, distribution-free alternative to Pearson, which is more appropriate for composite scores derived from ordinal Likert items (Boateng et al. 2018).

We visually inspect factor correlations as a heat-map and hide perfect correlations of a factor with itself to avoid visual anchoring on values of 1.0. High inter-factor correlations (above 0.7) strongly suggests the factors may be measuring the same underlying construct and should be considered for merging while medium/high correlations (0.5-0.7) indicate scrutinising whether the factors are conceptually distinguishable.

factors <- cols_tpsai %>%
  split(sub("_.*", "", .)) %>%
  imap_dfc(~ efa_data %>%
    transmute(!!.y := rowMeans(pick(all_of(.x)), na.rm = TRUE))
  )

cor(factors, method = "spearman", use = "pairwise.complete.obs") %>%
  as.data.frame() %>% 
  rownames_to_column("F1") %>%
  pivot_longer(-F1, names_to = "F2", values_to = "r") %>% 
  filter(F1 > F2)  %>% 
  ggplot(aes(x = F1, y = F2, fill = r)) +
  geom_tile() +
  geom_text(aes(label = round(r, 2))) +
  scale_fill_gradient2(
    low = "steelblue", mid = "white", high = "firebrick",
    midpoint = 0, limits = c(-1, 1)
  ) +
  labs(x = "", y = "", fill = "") +
  theme_minimal(base_size = 12) +
  theme(panel.grid = element_blank()) +
  scale_y_discrete(limits = rev)

Cronbach Alpha

Cronbach’s alpha is a measure of internal consistency reliability, quantifying the degree to which items within a factor covary and can therefore be considered indicators of the same underlying construct. Conceptually, it reflects the average inter-item correlation within a factor, weighted by the number of items. Importantly, alpha is not a measure of unidimensionality: a factor can have acceptable alpha while still containing sub-dimensions, which is why we complement it with item-level analyses below.

We then assess the internal consistency of each factor using ordinal Cronbach’s alpha, computed from the polychoric correlation matrix to account for the ordinal nature of our items.We interpret α ≥ .90 as excellent, .80–.89 as good, .70–.79 as acceptable, .60–.69 as questionable, and below .60 as poor.

cronbach_data <- cols_tpsai %>%
  split(sub("_.*", "", .)) %>%
  imap_dfr(~ {
    items <- efa_data %>% select(all_of(.x))
    rho   <- polychoric(items)$rho
    a     <- alpha(rho, n.obs = nrow(items))
    
    tibble(
      factor       = .y,
      alpha        = a$total$raw_alpha,
      ci_lower     = a$total$raw_alpha - 1.96 * a$total$ase,
      ci_upper     = a$total$raw_alpha + 1.96 * a$total$ase,
      n_items      = length(.x)
    )
  }) %>%
  mutate(
    interpretation = case_when(
      alpha >= .90 ~ "Excellent",
      alpha >= .80 ~ "Good",
      alpha >= .70 ~ "Acceptable",
      alpha >= .60 ~ "Questionable",
      TRUE         ~ "Poor"
    )
  )

Some items ( CRIT_02 ) were negatively correlated with the first principal component and 
probably should be reversed.  
To do this, run the function again with the 'check.keys=TRUE' option

Factors with questionable alpha are flagged for item-level inspection, as targeted removal of weak items may bring consistency to an acceptable level. Factors with poor alpha (< .60) are considered insufficiently coherent and candidates for removal entirely. At the upper end, factors with α > .95 are inspected for item redundancy, as excessively high consistency can indicate that items are too similar to justify their inclusion as distinct indicators.

cronbach_data %>%
  knitr::kable(
    digits  = 2,
    col.names = c("Factor", "α", "95% CI Lower", "95% CI Upper", "Items", "Interpretation")
  )

Factor	α	95% CI Lower	95% CI Upper	Items	Interpretation
COMP	0.85	0.82	0.88	5	Good
CORP	0.80	0.76	0.84	4	Good
CRIT	0.22	0.07	0.37	5	Poor
EMOT	0.77	0.73	0.82	5	Acceptable
HUMN	0.72	0.67	0.78	5	Acceptable
PERS	0.86	0.83	0.88	5	Good
RISK	0.78	0.74	0.83	5	Acceptable
SAFE	0.81	0.78	0.85	5	Good

Item-factor correlation

For each item, we compute the adjusted item-factor correlation, the correlation between that item and the sum of all other items in the same factor, excluding itself. Items with an adjusted item-total correlation below .30 contribute little to what the factor is seeking to measure (Boateng et al. 2018). Items with negative item-total correlations are of particular concern, as they suggest the item is measuring something in the opposite direction to the rest of the factor.

poly_r_drop <- function(items_df) {
  rho <- polychoric(items_df)$rho
  n   <- ncol(rho)
  
  map_dfr(seq_len(n), function(i) {
    others  <- setdiff(seq_len(n), i)
    cov_ir  <- sum(rho[i, others])
    var_rest <- sum(rho[others, others])
    tibble(
      item    = colnames(rho)[i],
      it_corr = cov_ir / sqrt(var_rest)
    )
  })
}

cols_tpsai %>%
  split(sub("_.*", "", .)) %>%
  map_dfr(~ {
    efa_data %>%
      select(all_of(.x)) %>%
      poly_r_drop()
  }, .id = "factor") %>%
  mutate(item_num = sub(".*_", "", item)) %>%
  ggplot(aes(x = factor, y = item_num, fill = it_corr)) +
  geom_tile() +
  geom_text(aes(label = round(it_corr, 2))) +
    scale_fill_gradientn(
      colours = c("#c94040", "#f0a875", "#efefef", "#82c882"),
      values  = scales::rescale(c(-1, 0, 0.3, 1)),
      limits  = c(-1, 1)
  ) +
  labs(x = "", y = "Item No.", fill = "r") +
  theme_minimal(base_size = 12) +
  theme(panel.grid = element_blank()) + 
  scale_x_discrete(position = "top") + 
  scale_y_discrete(limits = rev)

Cross-factor correlations

Item-factor correlations

We repeat the above analysis of item-factor correlations but this time looking at items correlate with factors that they are not assigned to it.

rho_full <- efa_data %>% 
  select(all_of(cols_tpsai)) %>% 
  polychoric() %>% 
  .$rho

factors <- split(cols_tpsai, sub("_.*", "", cols_tpsai))

imap_dfr(factors, function(focal_cols, focal_name) {
  imap_dfr(factors[names(factors) != focal_name], function(other_cols, other_name) {
    map_dfr(focal_cols, function(col) {
      idx_item  <- which(rownames(rho_full) == col)
      idx_other <- which(rownames(rho_full) %in% other_cols)
      cov_ir    <- sum(rho_full[idx_item, idx_other])
      var_rest  <- sum(rho_full[idx_other, idx_other])
      tibble(item = col, other_factor = other_name, it_corr = cov_ir / sqrt(var_rest))
    })
  })
}) %>% 
  mutate(item = factor(item, levels = rev(cols_tpsai))) %>%
  ggplot(aes(x = other_factor, y = item, fill = it_corr)) +
  geom_tile() +
  geom_text(aes(label = round(it_corr, 2)), size = 3) +
  scale_fill_gradientn(
    colours = c("#e8883a", "#efefef", "#c94040"),
    values  = scales::rescale(c(-1, 0, 1)),
    limits  = c(-1, 1)
  ) +
  labs(x = "", y = "", fill = "r") +
  theme_minimal(base_size = 12) +
  theme(panel.grid = element_blank()) +
  scale_x_discrete(position = "top")

Inter-item correlations

Lastly, we examine pairwise inter-item correlations. Items correlating above .7 suggest they may be measuring the same underlying concept. Where such pairs are identified, we consider whether the affected items should be merged, reassigned to a different factor, or removed.

ii_cor_threshold <- 0.7

efa_data %>%
  select(all_of(cols_tpsai)) %>%
  polychoric() %>%.$rho %>%
  as.data.frame() %>%
  rownames_to_column("item1") %>%
  pivot_longer(-item1, names_to = "item2", values_to = "r") %>%
  filter(r > -1 & r < 1) %>%
  filter(item1 < item2) %>% 
  filter(sub("_.*", "", item1) != sub("_.*", "", item2)) %>% # within-factor filter
  filter(abs(r) > ii_cor_threshold) %>%
  knitr::kable(digits = 2)

item1	item2	r

References

Boateng, Godfred O., Torsten B. Neilands, Edward A. Frongillo, Hugo R. Melgar-Quiñonez, and Sera L. Young. 2018. “Best Practices for Developing and Validating Scales for Health, Social, and Behavioral Research: A Primer.” Frontiers in Public Health 6 (June). https://doi.org/10.3389/fpubh.2018.00149.

Buchanan, Erin M., Kathrene D. Valentine, and Stefan E. Schulenberg. 2014. “Exploratory and Confirmatory Factor Analysis: Developing the Purpose in Life Testshort Form.” SAGE Research Methods Cases, January. https://doi.org/10.1007/S10902-008-9124-3.