Exploratory Factor Analysis Tables in R

The Factor Analysis

For this example, I am going to use Allison Horst’s palmerpenguins package for the data and produce a simple and meaningless factor structure. This is for demonstration purposes only

# load key packages
library(palmerpenguins)
library(tidyverse)
library(psych)
library(flextable)

# load data
data <- penguins %>%
  select(bill_length_mm:body_mass_g)

# factor analysis
efa <- fa(data, nfactors=4, rotate="oblimin")

Saving the factor analysis as efa creates a list with a number of useful elements:

names(efa)

##  [1] "residual"      "dof"           "chi"           "nh"           
##  [5] "rms"           "EPVAL"         "crms"          "EBIC"         
##  [9] "ESABIC"        "fit"           "fit.off"       "sd"           
## [13] "factors"       "complexity"    "n.obs"         "objective"    
## [17] "criteria"      "STATISTIC"     "PVAL"          "Call"         
## [21] "null.model"    "null.dof"      "null.chisq"    "TLI"          
## [25] "F0"            "r.scores"      "R2"            "valid"        
## [29] "score.cor"     "weights"       "rotation"      "communality"  
## [33] "communalities" "uniquenesses"  "values"        "e.values"     
## [37] "loadings"      "model"         "fm"            "rot.mat"      
## [41] "Phi"           "Structure"     "method"        "scores"       
## [45] "R2.scores"     "r"             "np.obs"        "fn"           
## [49] "Vaccounted"

For our purposes, we are going to focus on Vaccounted and loadings. Vaccounted is comprised of a table of eigenvalues and their respective variances accounted for:

efa$Vaccounted

##                             MR1       MR2        MR3          MR4
## SS loadings           2.1807053 0.5820772 0.11719588 4.667710e-17
## Proportion Var        0.5451763 0.1455193 0.02929897 1.166928e-17
## Cumulative Var        0.5451763 0.6906956 0.71999460 7.199946e-01
## Proportion Explained  0.7571950 0.2021117 0.04069332 1.620745e-17
## Cumulative Proportion 0.7571950 0.9593067 1.00000000 1.000000e+00

For this table, we are simply going to convert it to a data frame and then make it a little bit fancier with flextable. I have been including the flex function in a lot of my assignment files lately. flex is a function that takes a data frame and uses flextable to style it in a Word-friendly, APA-like format:

# table set-up ----

flex <- function(data, title=NULL) {
  # this grabs the data and converts it to a flextbale
  flextable(data) %>%
  # this makes the table fill the page width
  set_table_properties(layout = "autofit", width = 1) %>%
  # font size
  fontsize(size=10, part="all") %>%
    #this adds a ttitlecreates an automatic table number
      set_caption(title, 
                  autonum = officer::run_autonum(seq_id = "tab", 
                                                 pre_label = "Table ", 
                                                 post_label = "\n", 
                                                 bkm = "anytable")) %>%
  # font type
  font(fontname="Times New Roman", part="all")
}

I use the flex function to print tables that are automatically numbered in word. You will need to add the following to your YAML to make that work:

output: bookdown::word_document2

Returning to our table, the following code works easily:

efa[["Vaccounted"]] %>%
  as.data.frame() %>%
  #select(1:5) %>% Use this if you have many factors and only want to show a certain number
  rownames_to_column("Property") %>%
    mutate(across(where(is.numeric), round, 3)) %>%
    flex("Eigenvalues and Variance Explained for Rotated Factor Solution")

Simple enough, right?

Now, on to a more slightly more complicated table. The loadings are stored in the loadings element, but we also might want to communalities, uniqueness, and complexity score for each item too. Those are stored in different elements. In addition, we will want to sort the factor analysis so that it is easy to read. We’ll need to use fa.sort for that. We also want to show factor loadings and suppress anything under .32 (a standard cut-off value). Finally, we will want to make this a simple function call so that we can use it multiple times with easy - factor analysis is an iterative process, after all. Here is the fa_table function:

fa_table <- function(x, cut) {
  #get sorted loadings
  loadings <- fa.sort(x)$loadings %>% round(3)
  #supress loadings
  loadings[loadings < cut] <- ""
  #get additional info
  add_info <- cbind(x$communalities, 
                    x$uniquenesses,
                    x$complexity) %>%
    # make it a data frame
    as.data.frame() %>%
    # column names
    rename("Communality" = V1,
           "Uniqueness" = V2,
           "Complexity" = V3) %>%
    #get the item names from the vector
    rownames_to_column("item")
  #build table
  loadings %>%
    unclass() %>%
    as.data.frame() %>%
    rownames_to_column("item") %>%
    left_join(add_info) %>%
    mutate(across(where(is.numeric), round, 3))
}

And here is the result:

fa_table(efa, .32)

## Joining, by = "item"

##                item   MR1   MR2 MR3 MR4 Communality Uniqueness Complexity
## 1       body_mass_g 0.932                     0.799      0.201      1.030
## 2 flipper_length_mm 0.791                     0.995      0.005      1.252
## 3    bill_length_mm 0.718                     0.584      0.416      1.311
## 4     bill_depth_mm       0.606               0.502      0.498      1.126

And here is the table:

fa_table(efa, .32) %>%
  flex("A Pretty Factor Analysis Table")

Any factor loadings under .32 were suppressed, meaning they are an empty character string that shows as blanks. Had this been a long table with many more factors, the default output in R would have been unwieldy to view and difficult to print. Using the fa_table function makes it easier.