Factor Analysis / PCA

A personality survey might ask 50 questions, but the answers are not 50 independent things. They cluster, because the questions tap a handful of underlying traits. Factor analysis and Principal Component Analysis (PCA) find those few dimensions among many correlated variables, so that 50 columns can be summarized by perhaps 5.

The core idea: directions of variance

PCA looks for the directions in which the data varies the most. The first principal component is the single line that captures as much of the spread as possible. The second is the best remaining direction perpendicular to the first, and so on. When variables are highly correlated, almost all the variance lines up along one direction, so one component can stand in for many.

Watch the last readout. Two variables span two dimensions, but when they are strongly correlated the cloud is close to a line, so it is nearly one-dimensional. PCA measures this: if PC1 explains 95% of the variance, one number per point describes the data almost as well as two. With 50 variables, this is dimensionality reduction, and the same method reappears in machine learning.

PCA vs. factor analysis

The model behind factor analysis

Factor analysis is built on a measurement model, an equation for how an observed answer is produced. For person n answering item i:

y_i,n = μ_i + β_i ξ_n + ε_i,n

Read it left to right. An answer starts from the item's own mean (μ_i, mu). To that is added how strongly the item responds to the trait (β_i, beta, its loading), multiplied by how much of the trait this person has (ξ_n, xi, their factor score). Last comes whatever that answer had of its own (ε_i,n, epsilon, the item's error).

ξ is never observed, so it has no units of its own. Something has to fix its scale, and the usual convention is to standardize the factor to mean 0 and variance 1. Every factor loading is then on the same scale, so software can report loadings without asking what units the factor is measured in. Loading tables are normally computed on standardized items too, and then each β is the correlation between that item and the factor.

Two identities follow, and the rest of this lesson builds on them. The first is about two items. The errors are assumed uncorrelated with the factor and with each other, so items i and j are related only through ξ. Their covariance (σ, sigma, which on standardized items is their correlation) is therefore the product of their two loadings:

σ_ij = β_i β_j

The second is about one item. Its variance splits into the part the factor explains and the part that belongs to the item alone, written ψ_ii (psi):

σ_ii = β_i² + ψ_ii

Each part of that split has a name. β_i² is the item's communality, written h²: the share of its variance the factor accounts for. ψ_ii is its uniqueness, partly content specific to the item and partly measurement error. The share of an item's variance that the factor accounts for is also what reliability means for a single indicator. So with one factor and standardized items, an item's communality is its reliability. The two coincide because σ_ii = 1 after standardizing. On raw, unstandardized items, reliability is β_i² ÷ σ_ii, and the two differ.

For example, three items loading .80, .70 and .50 on one factor imply correlations of .56, .40 and .35 between them, communalities of .64, .49 and .25, and uniquenesses of .36, .51 and .75. Each communality and its uniqueness add to 1, because standardizing gave each item exactly one unit of variance.

Both identities are statements about covariances, and a correlation is a covariance divided by the two standard deviations. So the model says nothing about the item means: μ_i drops out as soon as you center the items. A factor analysis therefore runs on the p × p matrix of correlations rather than on the N × p raw data. That is why you can reproduce a published factor analysis from the article's correlation table, and why SPSS will run one from a matrix typed in by hand.

A factor solution is judged by its residual grid. At the fitted loadings every residual is under .025 in absolute value (the largest is −.024), so four loadings reproduce the six observed correlations almost exactly. Drag β₁ down to .40 and three residuals in q1's row jump past .20 at once, because that one loading appears in three of the six products. This is what "the model reproduces the correlation matrix" means. Extraction searches for the loadings that bring this grid as close to zero as possible, and the goodness-of-fit test further down asks whether what is left is more than sampling noise.

The root mean square residual in the readout summarizes that grid in one number: the square root of the average squared off-diagonal residual. Below about .05, a solution is generally considered to have reproduced the matrix, and the fitted values here give .013. Read the grid as well, though. One large residual among tiny ones means a pair of items shares something the factor does not capture, and the cause is usually the wording of those items rather than the model.

Is the matrix worth factoring at all?

Before extracting anything, check whether your variables are correlated enough to have shared structure. Software offers two standard checks.

KMO (the Kaiser-Meyer-Olkin measure of sampling adequacy) compares the ordinary correlations between your variables with their partial correlations, the ones left after every other variable has been held constant. When the variables share common factors, holding the others constant removes most of each correlation, so the partial correlations are small and KMO climbs toward 1. Variables that correlate only in isolated pairs keep large partial correlations and push it down. Kaiser's labels run from marvelous above .90 through middling in the .70s to unacceptable below .50, and .60 is the usual working floor. A clean two-factor dataset scores around .91, and pure noise about .50.

Bartlett's test of sphericity asks a much weaker question: could your correlation matrix have come from variables that are all mutually uncorrelated? Treat a significant result as a formality rather than an achievement. With 20 variables and 100 people, true correlations of just .10 already give p = .00005, and at 300 people the same tiny correlations give p below 10⁻¹⁶. A non-significant Bartlett test means stop. A significant one means almost nothing on its own, so rely on KMO and the scree plot.

Cronbach's α (alpha) is usually run in the same session, but it answers a different question. KMO and Bartlett ask whether the matrix has any shared structure. α assumes the items belong together and asks how consistently they behave as one scale:

α = (k / (k − 1)) × (1 − Σσ_ii ÷ σ²_total)

k is the number of items, σ_ii each item's variance and σ²_total the variance of the total score. The ratio inside the bracket is the share of total-score variance that comes from the items' separate variances. The rest comes from their covariances, and that shared part is what α reports. A floor of .70 is conventional and .80 is comfortable for research use. Above about .95 is a warning, because it usually means several items are near-duplicates of each other. The eight wellbeing items score .90.

Two limitations of α matter when a reviewer asks about it.

The first is that α rises with the number of items, whatever their quality. Hold the average correlation between items at a mediocre .20 and keep adding items:

The items did not improve along that row, yet thirty items correlating .20 with each other reach a respectable .88. The Spearman-Brown formula states the same logic directly: averaging more indicators cancels more item-specific noise. Report α beside k, and do not take a high α on a long scale as evidence that the items are good.

The second is that α is about internal consistency, not about unidimensionality. It does not tell you how many factors there are. A questionnaire built from two clearly separate factors, with items that correlate reasonably within each block, can still return an α in the .80s. So a high α is no evidence that a scale measures only one factor. Use the factor analysis to establish the scale's structure, and α to report how consistently the established scale behaves.

Choosing an extraction method

Extraction turns a correlation matrix into loadings. Software offers several methods, and two of them account for nearly everything you will read.

Principal axis factoring applies the identities above directly. It puts an estimate of each item's communality on the diagonal of the correlation matrix in place of the 1s and solves for the loadings. Then it uses those loadings to re-estimate the communalities, puts the new values on the diagonal, and repeats until they stop changing. It converges on the loadings that make the off-diagonal residuals as small as possible, the same thing you were doing by hand in the widget above. Principal axis factoring assumes nothing about the shape of the data, and that is the reason to choose it: Likert answers are ordinal, often skewed, and rarely normal.

Maximum likelihood asks a different question: if the variables are multivariate normal, which loadings make the observed data most likely? The normality assumption is a real cost. In return you get a likelihood, and with it a formal goodness-of-fit test of whether the number of factors you asked for is enough. Principal axis factoring has no such test, so the worked χ² further down this page comes from maximum likelihood.

Choose principal axis factoring when normality is doubtful and you want a defensible structure. Choose maximum likelihood when the items are plausibly continuous and roughly normal and you want to test a specific number of factors. In SPSS both are in one dropdown under Extraction…. The default there is neither of them, so always set the method yourself.

How many factors to keep

An eigenvalue is usually defined as "the variance a factor accounts for". In terms of loadings, a factor's eigenvalue is the sum of its squared loadings down the column. Each β_i² is the share of item i that the factor explains, so adding them across the items counts how many whole items' worth of variance the factor accounts for. In the one-factor solution for the eight wellbeing items, the squared loadings run from .444 to .575 and total 4.32. So one factor accounts for as much variance as about four and a third of the eight items, 54% of the eight units of variance there are to explain.

This explains the Kaiser criterion (keep factors with an eigenvalue above 1). Every item is standardized to one unit of variance, so 1 is the variance of a single item. A factor with an eigenvalue below 1 explains less variance than one of your own questions does, and you would be better off keeping the question.

The scree plot draws those eigenvalues in descending order, and you look for the "elbow" where the line flattens. A third habit is to keep enough factors to reach a target share of the variance. The target is often 80 or 90% in a components analysis and lower in a factor analysis, because a factor analysis does not try to explain the unique variance.

Parallel analysis is the retention rule most methodologists now recommend. Simulate many random datasets with your number of variables and your number of people, average their eigenvalues, and keep only the factors whose real eigenvalue is larger than the matching random one. The comparison is needed because 1 is a low bar. With 12 variables and 100 respondents of pure noise, the largest eigenvalue averages 1.61, and on average 5.6 of the 12 eigenvalues exceed 1.0. The Kaiser criterion would keep factors built entirely from sampling variation.

Expect the rules to disagree. Scree and Kaiser often give different counts on the same matrix, and parallel analysis often gives a count between the two. No further statistic settles it. Decide by interpretability and theory: fit each candidate number of factors, rotate each solution, and keep the one whose factors you can name and defend to someone who knows the construct. A two-factor solution you can explain is better than a three-factor solution whose third factor is two items you cannot interpret.

Does the model fit?

The fit test that maximum likelihood provides needs careful reading, because its null hypothesis runs the opposite way to almost every other one on this site. The null says that m factors are enough: the model-implied correlations are the population's own, and every residual you see is sampling noise. So a small p is bad news here, and you are hoping for a result you cannot reject.

The test's degrees of freedom come from counting. A correlation matrix over p items holds p(p + 1) ÷ 2 distinct numbers. The model estimates pm loadings and p uniquenesses, but m(m − 1) ÷ 2 of the loadings are redundant, because rotation can change them without changing the fit. Subtract and collect terms:

df = ((p − m)² − (p + m)) ÷ 2

For the eight wellbeing items, the df fall quickly as m rises:

An eight-item correlation matrix holds 36 numbers, and one factor uses 16 parameters, which leaves 20 df for this page's one-factor fit: χ²(20) = 20.14, p = .45. One factor reproduces those correlations as well as sampling noise would allow. Three factors leave only seven df. With enough factors the model has as many parameters as the matrix has numbers, and then it fits perfectly and tells you nothing.

Two cautions apply to the test. It exists only under maximum likelihood, so a principal-axis solution has no fit statistic. A write-up that reports one alongside "principal axis factoring" has taken it from a different analysis. And χ² grows with sample size, which this page's own data can show. Measured on 200 people instead of 75, the same eight-item discrepancy would return χ²(20) = 56.19 and p < .001, and on 1,000 people χ²(20) = 286.91. Nothing about the fit changed. It is the same problem as a normality test rejecting a harmless departure in a large sample, with the same answer: read the test alongside the residual grid. A significant χ² with every residual under .05 is a large-sample artifact. A significant χ² with three residuals past .15 means the model is missing something.

Rotation: the same model, made readable

The solution that comes out of extraction is usually impossible to interpret. The first factor is a general factor that every item loads on, because extraction puts as much variance as it can into the first factor. The second factor describes whatever separates the items, which makes it bipolar: half the items load positively and half negatively. That describes the data accurately but does not match the constructs. Rotation turns it into a solution that does.

Rotation is allowed because a factor solution is not unique, and the algebra takes two lines. Write Λ (a capital lambda) for the matrix of loadings. Everything the model says about the data is in Λ Λ^T, the grid of implied correlations from earlier in this lesson. Now take any orthogonal matrix T, one whose own transpose is its inverse, and rotate the loadings with it:

(ΛT)(ΛT)^T = Λ T T^T Λ^T = Λ Λ^T

T and its inverse cancel in the middle. Every implied correlation, every communality, every residual and the total variance explained come out identical, while the loading table can look completely different. The same cancellation works for any invertible T, not only an orthogonal one, if the factors are allowed to correlate. Put Λ* = ΛT⁻¹ and let Φ (phi, the matrix of correlations between the factors) be TT^T, and Λ* Φ Λ*^T is Λ Λ^T again. This is the whole justification for rotation, and it also sets its limit. Every rotation fits the data equally well, so rotation cannot improve a solution. It can only make it easier to interpret.

That leaves a choice about Φ, and the choice is a claim about your constructs. Orthogonal rotation, of which varimax is the most widely used, holds Φ at the identity: the factors are forced to be uncorrelated. Oblique rotation, whose common forms are promax and direct oblimin, lets them correlate and then reports how much. In psychology an oblique rotation is the usual recommendation. A theory about two related constructs does not claim they are independent, and an oblique rotation estimates how strongly they correlate. Run oblimin first and read Φ. If the factor correlation comes out near zero, varimax is defensible after all, and you have evidence for that choice.

At 0° the solution is the one extraction produced, as described above. Every item loads between .53 and .76 on the first factor, the second factor is bipolar, and only one of the eight items loads on a single factor by the ±.32 rule. Turn the dial to −44° and seven of the eight do. Nothing about the model changed. The last readout is the largest change in any of the 28 implied correlations, and it stays at .000 at every angle in both modes, because the cancellation above is exact.

The first readout is the criterion varimax maximizes: the variance of the squared loadings. It is large when each column holds a few big values and many near-zeros, and small when every loading is middling. Each item's row is scaled to unit length before the variance is taken, so a well-explained item does not count for more than a poorly explained one. That scaling is Kaiser normalization, and it is on by default in SPSS. The software finds the peak you can search for by hand with the dial, and the button goes to the same angle, −44° for these items. Signs carry no meaning: a factor can be multiplied by −1 without changing anything, and software sometimes reports a solution with all the signs flipped relative to the paper you are replicating.

Pattern or structure: which table do I read?

After an orthogonal rotation there is one rotated factor matrix and no ambiguity. After an oblique one, SPSS prints two, and they are not the same numbers.

The pattern matrix holds each factor's own contribution to an item with the other factors held constant, so a pattern coefficient is a regression weight, not a correlation. The structure matrix holds plain correlations between items and factors. When the factors are correlated, those correlations are inflated: an item that belongs to factor 2 will also correlate with factor 1, only because the two factors are correlated. Interpret the pattern matrix. Reading the structure matrix instead is a common mistake, and the widget above shows why. Switch it to oblique and press the button: it stops at −32°, with the factors correlating .40. In the pattern columns, seven items load on exactly one factor. In the structure columns, on the same solution and the same data, only four do. Items i1, i5 and i7 each gain a second loading past .32, although their pattern coefficients on that factor are .002, .048 and .076. Those three cross-loadings come entirely from the correlation between the factors.

One consequence often catches students out. Under an orthogonal rotation, an item's communality is the sum of its squared loadings. Under an oblique one it is not, because the factors overlap and the overlap would be counted twice. The communality column does not change, since rotation cannot change it, but you can no longer get it by squaring and adding the pattern coefficients.

Reading the loading table

Two thresholds do most of the work. A loading of .32 is the usual minimum for an item to count toward a factor: .32² = .1024, so the factor explains about a tenth of that item's variance. A loading of .50 squares to .25, a quarter of the item's variance, and is the conventional mark of a strong indicator. Some fields use .40 as the minimum instead. State which threshold you are using before you apply it, not after seeing which one tidies your table.

Reverse-code negatively worded items before the analysis, not after. On a 1 to 5 agreement scale the reversed score is 6 − x, so 1 becomes 5, 2 becomes 4 and 3 stays the same. On a scale running 1 to k, the general form is (k + 1) − x. If you skip this step, the item correlates negatively with everything else in its block. That does more than flip a sign in the output: extraction can produce a second factor made entirely of the reversed items. The wellbeing data on this page has two such items, q3 and q6, and both are reverse-coded before anything else. Recode into a new variable rather than over the original, so the raw data survive if you make a mistake.

A cross-loading has two conditions, and people often forget the second. An item cross-loads when it loads above the threshold on two factors and those two loadings are within about .20 of each other. Item i4 in the widget is an example: after varimax it reads .589 and .485, both past .32 and .104 apart, so it does not clearly belong to either factor. An item reading .62 and .35 is different. It is over the threshold twice, but it is a factor-1 item with a smaller secondary loading, and the write-up can say so.

When an item does cross-load, work through these steps in order.

1. Re-specify the number of factors. A cross-loading is often a sign that the number of factors is wrong by one, and the item loads cleanly once the right number is fitted.
2. Try the other rotation. An item between two correlated factors may load cleanly under an oblique rotation, where an orthogonal rotation forced it to split between them.
3. Drop the item, and say in the write-up that you did and why. This is a legitimate option, but it comes third, because the information in a dropped item is lost.
4. Check what is left. Each factor needs at least three indicators to be stable. A factor with two is fragile, and one with a single item is not a factor at all. If dropping the item leaves a factor with two, the problem is the questionnaire, not the item.

Then name the factors. No statistic helps with this. Read the items loading highest on each factor and write a label that a colleague could predict from the items alone. If you can name a factor only by listing its loadings, you have not interpreted it.

How much of an item does the model explain?

With several factors, an item's communality is the sum of its squared loadings across the factors you kept, in an unrotated or orthogonally rotated solution. The uniqueness is still what is left. An item loading .70 on a single factor has h² = .49, so the model explains just under half of its variance.

Fit a one-factor model to the eight wellbeing items in wellbeing-survey.csv (after reverse-coding q3 and q6) and the loadings are between .67 and .76. The communalities run from .44 to .58, with a mean of .54. Those same eight items have a Cronbach's α of .90. Both numbers are correct, and together they make a useful point. Any single item is about half noise, but their average is not, because the item-specific halves are unrelated to each other and largely cancel. That is why a scale of several items is more reliable than one good question, and it is reliability seen from the item's side.

Communality is also where the arithmetic of the two methods differs. PCA analyzes the correlation matrix with 1s down the diagonal, so it treats every variable as fully explainable and repackages all of the variance, error included. Factor analysis puts estimated communalities on that diagonal, so it models only the variance the items share. With strong items the two solutions look almost identical. With weak ones, PCA reports larger loadings and more variance explained than the data can support.

Read the communality column before dropping an item for a low loading, because an item can load modestly on each of three correlated factors and still be well explained. A communality near .20, on the other hand, says something about the item rather than the model: whatever that item measures, the factors you kept do not contain it. If the items came from a questionnaire, the fix is usually in how the item was worded, not in the extraction settings.

Factor scores

Often you want a number per person rather than per item: a wellbeing score to correlate with something else, or to use as a predictor in a regression. That number is a factor score, and the precise word for getting one is predict, not compute or estimate. You compute a quantity that is fully determined by the data, like a mean. You estimate a parameter, a fixed unknown constant of the population, like a loading. A person's ξ_n is neither. It is a real quantity that belongs to that person and was not observed. In that respect it is like a residual in a regression, which is also predicted rather than estimated.

Because ξ was never observed, there is no single right answer, and software offers several predictors. You will mostly meet two.

The regression predictor (SPSS calls it the regression method, and it is the default) chooses the weights that minimize the mean squared difference between the predicted scores and the true ones. In that sense it is the most accurate predictor available. It also takes the correlations between factors into account, which matters after an oblique rotation. The cost is that its scores are not specific to one factor. A score for factor 1 also correlates somewhat with factor 2, so two factor scores from a correlated pair are more correlated than the factors themselves.

The Bartlett predictor makes the opposite trade. Its scores are unbiased: across repeated samples they center on the true factor score rather than shrinking toward the mean. And each score correlates with its own factor and not with the others. Use it when the scores themselves are the outcome of interest and you need to say that this column measures this factor and nothing else.

After an orthogonal rotation either is defensible, and the choice depends on whether you care more about accuracy or unbiasedness. After an oblique rotation the regression predictor is the usual choice, because it takes the factor correlations into account. Report which one you used: the two give different scores, and a reader cannot tell which from the numbers.

Whichever you choose, compare it with the simple alternative. A plain sum or average of the items usually correlates above .95 with the predicted scores, works in a new sample without refitting, and is what most published scales report. On the eight wellbeing items the two agree at r = .999, as they normally do when the loadings are similar. Predicted scores are the better choice when the loadings are clearly uneven or the factors correlate. Otherwise they are not better than adding the items up.

The order you do it in

The analysis runs in six steps, in the order a practical exam usually asks for them.

1. Explore the data first. Reverse-code the negatively worded items, look at the correlation matrix, and run KMO and Bartlett's test to check there is shared structure to find at all.
2. Decide how many factors. Use the scree plot, the Kaiser criterion and parallel analysis, with theory to break ties, plus the fit test if you extracted by maximum likelihood.
3. Estimate the loadings, naming the extraction method and why: principal axis factoring when normality is doubtful, maximum likelihood when it is not and you want the fit test.
4. Rotate. Orthogonal if you are claiming independent constructs, oblique if you are not, and read Φ to find out which claim the data supports.
5. Interpret and label. Apply your stated loading threshold, deal with the cross-loadings, check each factor still has three indicators, and give the factors names.
6. Predict factor scores, if you need a number per person for what comes next.

Steps 2 and 3 are not as sequential as the list makes them look. Fitting two candidate numbers of factors and comparing the rotated solutions is normal practice, provided you report that you did it.

Problem 81 of the practice problems gives three loadings and the correlations the three items produced. It asks for the implied matrix, the communalities and the residuals, then what happens when you refit one loading. Problem 82 works through α's dependence on the number of items: a six-item scale, the same six items written twice over, and a reliability coefficient that rises without any item improving. Problem 83 takes a three-factor extraction apart with eigenvalues and the fit test, including a factor that explains only a sixth of one item's worth of variance. Problem 84 is a rotated oblique solution to label, with a structure matrix that makes every item look like a cross-loading, though only one item is.