skip to main content

A Factor Analysis For Dummies in R

published icon  |  category icon education software

tags icon phd R factor analysis

Suppose you want to develop a survey, like we recently did for creativity in higher computing education. You come up with a bunch of questions (Time seemed to fly while programming!, I adapted and merged ideas from others, etc), recruit victims who are willing to answer them using a scale (Totally agreed, Agreed, No way, etc), and then what? Chances are you picked certain questions from certain subthemes because you suspect a pattern will emerge. How to validate that?

That’s where “factor analysis” comes in: it’s a statistical method—meaning it can analyze quantitative data, like answers from the above scale—to detect variability and correlations among the questions. My research is very much qualitative (looking for patterns using interviews and discussions), so I had no idea on how to handle that. This blog post servers as a summary for myself and can perhaps entertain others.

Step 1: Is the data up for it?

First things first: it’s possible that your questions are so bad there’s never going to emerge any pattern. In that case, the “factorability” of the survey is reduced to zero.

There’s a statistic for that (That sentence comes in handy later!), it’s called the Kaiser-Meyer-Olkin (KMO) sampling adequacy. Their original paper suggests a value of at least .60, but preferably higher than .80.

Let’s use a well-known sample data set for our exercise: the psychological Big Five personalty test developed in the eighties. Scientists devised five big personality types, exactly using this very method: openness (O), conscientiousness (C), extraversion (E), agreeableness (A), neuroticism (N). The psych package contains a sample dataset of 2800 responses that can be loaded into the variable bfi with data('bfi').

I’ll deliberately avoid using statistical formulas here and let R Suite do the talking:

library('psych') #import fa functions
data('bfi')
KMO(bfi)

The result:

Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = bfi)
Overall MSA =  0.84
MSA for each item = ...

All right, so far so good. What’s next? An overall item-total correlation (Rit) check to validate the reliability of the questions asked (the items). There are, of course, standardized rules for that, laid out by Ebel and Frisbie. Rit values below .30 are considered doubtful or even poor. How are we doing on that front? Let’s try out the alpha(bfi) command and inspect the results:

Reliability analysis   
Call: alpha(x = bfi)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
      0.21      0.51    0.66     0.036 1.1 0.019  4.6 0.53    0.017

 lower alpha upper     95% confidence boundaries
0.18 0.21 0.25 

 Reliability if an item is dropped:
          raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r med.r
A1             0.24      0.53    0.68     0.041 1.15    0.018 0.035 0.021
A2             0.19      0.49    0.64     0.034 0.95    0.019 0.034 0.017
A3             0.19      0.49    0.64     0.034 0.94    0.019 0.033 0.017
...

 Item statistics 
             n  raw.r  std.r  r.cor r.drop mean    sd
A1        2784 -0.028 0.0573 -0.057 -0.127  2.4  1.41
A2        2773  0.285 0.3736  0.356  0.210  4.8  1.17
A3        2774  0.267 0.3874  0.379  0.184  4.6  1.30
A4        2781  0.269 0.2938  0.239  0.174  4.7  1.48
A5        2784  0.265 0.3345  0.314  0.185  4.6  1.26
C1        2779  0.204 0.2723  0.215  0.125  4.5  1.24
C2        2776  0.229 0.3506  0.324  0.142  4.4  1.32
C3        2780  0.184 0.2487  0.181  0.099  4.3  1.29
C4        2774  0.040 0.1285  0.060 -0.055  2.6  1.38
C5        2784  0.083 0.1242  0.049 -0.027  3.3  1.63
E1        2777  0.032 0.0023 -0.107 -0.078  3.0  1.63
...

Lots of numbers, where to look? First, the raw_alpha is a total alpha value that should be as close to one as possible. A low alpha value (.21 is not good) predicts that some questions can be left out while retaining the core and purpose of the survey.

The r.cor column is what we’re interested in. Surprisingly, there are quite a few bad values here. For instance, A1 could be scrapped, C5 and E1 are also worthless. At this point it’s up to the researcher to decide whether or not to throw out questions. We did for a few, after carefully inspecting the question itself. Some participants might have misinterpreted the question. Of course, just throwing out all bad ones to up the statistical relevance is not the way to go.

You can re-match the questions to the codes in the R documentation: A1 is “I Am indifferent to the feelings of others.”, E1 is “I Don’t talk a lot.” Since we already know the 5 factors it’s a bit difficult to judge whether or not those questions are relevant to the personality types. Anyway, let’s move on and pretend we didn’t see the low alpha value.

Step 2: Confirm your suspicions

There are two ways to do a factor analysis: confirmatory or exploratory. With the first, you suspect certain items will belong together, and hope that the statistics will confirm that. With the latter, you have no idea what will emerge, and let the analysis tell you which ones should be grouped together.

Usually, researchers have a pretty good idea, especially if they designed the questions themselves. They’re just looking for confirmation—and a way to publish their results. Lets pretend we think O1, C1, E1, A1, N1 and O2, C2, … belong together, comprising four factors in total. We have to feed R our model ourselves, and then pass it to the cfa function from the lavaan project:

library('lavaan')

model <- '
 # measurement model
  cool  =~ O1 + C1 + E1 + A1 + N1
  mad =~ O2 + C2 + E2 + A2 + N2
  sad =~ O3 + C3 + E3 + A3 + N3
  glad =~ O4 + C4 + E4 + A4 + N4
' 

fit <- cfa(model, bfi)
summary(fit, fit.measures = TRUE)

The summary function spews out about every imaginable statistic you can think of. We’re looking for a few things here:

  • The Comparative Fit Index (CFI) that should ideally be >.80. Here, it’s .433.
  • The Root Mean Square Error of Approximation (RMSEA) and its P-value that should be below .05. Here, it’s .135.

All other stuff (variances, covariances, latent variables) are useless if the CFI is low, as is the case. That means your model is wrong, and the confirmatory factor analysis confirmed that your hunch is… well… just a hunch. What’s next, scrapping the data and having another go at collecting “better” ones? Not yet. It just means that our cool, mad, sad, glad factors are hogwash.

Step 3: Explore other possibilities

At this point, it’s a good idea to trade confirmatory for exploratory, although fiddling with the measurement model yourself is always an option—in my case, it only resulted in frustration. R still asks the researcher for an educated guess when it comes to the number of factors, though.

If you have no clue where to begin, a parallel analysis is a good idea, which comes with a bunch of optional parameters, including which factor method to use, which I won’t go into here:

eigenvalues <- fa.parallel(bfi)

This will take a while, but you can limit the result set by using bfi[1:100].

The result is one line of output (“Parallel analysis suggests that the number of factors = 7 and the number of components = 7”) and a graph that is called a “scree plot”:

Plotted eigenvalues using fa.parallel.

The x axis contains the number of potential factors, and the y axis the “eigenvalues”, which explains how much of the variance of the items a factor explains. An eigenvalue of one (the threshold) means that factor on the x axis explains more variance than a unique variable. In other words, If it goes below one, introducing new factors won’t be very beneficial, as the data of the questionnaire contributes little to that factor.

These graphs traditionally resemble a kind of elbow: after x factors, adding another one becomes futile. This is basically R telling us how many factors we should be aiming for. In this case, it could be seven, as suggested, but it could also be six or five. Remember: it’s still an educated guess.

Next, let’s try out that guess. Seven? Sure, why not.

factordata <- fa(bfi, 7, rotate="varimax")
> print(factordata$loadings, cutoff=0.3) 

The result:

Loadings:
          MR1    MR2    MR3    MR4    MR5    MR6    MR7   
A1                                    -0.538              
A2         0.391                       0.513              
A3         0.524                       0.434              
A4         0.302                       0.311              
A5         0.576                                          
C1                       0.544                            
C2                       0.670                            
C3                       0.542                            
C4                      -0.609                            
C5                      -0.537                            
E1        -0.504                              0.359       
E2        -0.597                              0.386       
E3         0.657                                          
E4         0.683                                          
E5         0.510         0.308                            
N1                0.796                                   
N2                0.794                                   
N3                0.723                                   
N4                0.551                       0.314       
N5                0.518                                   
O1         0.346               -0.393                     
O2                              0.517                     
O3         0.428               -0.498                     
O4                                            0.350       
O5                              0.591                     
gender                                 0.366              
education                                            0.389
age                                                  0.564

Whoops, what are those last few items? Demographic info can be removed but I was too lazy to do that. What are we looking at here? The aim for these “loadings” is to (1) have every item (question) loaded on a separate factor, and not multiple, and (2) have the questions grouped in the factors in a way that makes sense. MR1, MR2, etc can afterwards be named. Since we’ve kind of spoiled the ending, this is a bit ridiculous now: obviously, we’re looking to load all the A questions under a factor (“agreeableness”), all the Cs, etc.

Now is a good time to re-calculate that comparative fit index: .97. Wow, that’s great! Of course, since R basically suggested seven factors.

What does rotate=varimax mean? There are different kinds of matrix rotations available, grouped into orthogonal and oblique. As the analysis looks for the correlations to determine the best fit, it can rotate around an axis to find a better fit. In short, if you assume your factors are going to be correlated (which is mostly the case in psychology!), choose an oblique rotation. If you think they’re not, choose an orthogonal one.

In the above loading output, you can see that for instance E5 loads high on MR1 (.510) and MR3 (.308). That’s annoying. Is this question related to construct 1 or 2? Perhaps five factors is a better fit. Let’s try again with five and apply an oblimin rotation, because we suspect they’re correlated:

Loadings:
          MR2    MR1    MR3    MR5    MR4   
A1                             -0.514       
A2                              0.631       
A3                              0.581       
A4                              0.397       
A5                0.320         0.458       
C1                       0.547              
C2                       0.662              
C3                       0.563              
C4                      -0.623              
C5                      -0.560              
E1               -0.539                     
E2               -0.641                     
E3                0.490                     
E4                0.669                     
E5                0.407                     
N1         0.782                            
N2         0.762                            
N3         0.729                            
N4         0.502 -0.356                     
N5         0.531                            
O1                                     0.522
O2                                    -0.436
O3                                     0.616
O4                                     0.366
O5                                    -0.522

Bingo. Everything seems to be (almost) neatly grouped into five different factors—the ones we started with: openness (O), conscientiousness (C), extraversion (E), agreeableness (A), neuroticism (N). Of course, when scientists were originally analyzing this data, they didn’t know that yet, and had to come up with those cool labels. What about our CFI? .89—still great.

A few other functions that can be of use:

  • principal, calculating aforementioned statistics for an exploratory analysis
  • fa.diagram(bfi) visualizes the loaded factors

If you’re still keen on grinding away in R, you can verify the EFA findings by creating a revised model and doing a second CFA. The results should more or less match, including the high CFI value.

Step 4: Interpret with caution

What’s the point of all those statistics? The emerging factors are never 100% the expected factors. For instance, we expected to see the same themes as the ones we identified in a focus group, as the questions were based on interviewees' opinions, supplemented with academic literature. We instead found (1) less factors, (2) a unique group of sub-themes we didn’t foresee, and (3) discrepancies in our questions.

That means the process of coming up with questions, gathering data, and analyzing them using something like a factor analysis, is an incremental process. It should—in theory—be done multiple times using revisions of your questionnaire. In practice, gathering test subjects just once already involves a lot of trouble. I’ve found a lot of psychometric tests that never receive a published revision.

One more thing: to get your stuff really, really correct, you need a Step 0: Validate beforehand, or a “face validity” phase. There, you present your survey to a small group of experts/critics/potential candidates to iron out the biggest syntactic mistakes. That’s the only qualitative part of a survey design.

Now you know the basics of setting up and validating your own questionnaire. I had no idea what I was doing, and I still have no clue on some of the statistical calculations behind these functions, but at least now I know how to use them and interpret their results. Good luck!

I'm Wouter Groeneveld, a Brain Baker, and I love the smell of freshly baked thoughts (and bread) in the morning. I sometimes convince others to bake their brain (and bread) too.

If you found this article amusing and/or helpful, you can buy me a coffee - although I'm more of a tea fan myself. I also like to hear your feedback via Mastodon or e-mail. Thanks!