Principal components and factor analysis for marketers

TL;DR

Difficulty: Advanced

  • Twenty four correlated survey questions are usually four or five underlying things measured many times. PCA and factor analysis find those things.
  • PCA turns the data so the variance sits at the front and assumes nothing. Factor analysis is a model in which hidden constructs cause the answers and each item keeps a unique part.
  • Use factor analysis for surveys that measure constructs, PCA to compress behaviour before clustering or regression. Standardise first, log skewed variables.
  • Decide the count with the scree plot, Kaiser's rule and parallel analysis together, and let interpretability overrule them. Rotation changes readability, never fit.
  • Scores only pay off once they are tied to an outcome: in the example, fair pricing of extras beats a room refurbishment.

Most surveys and customer tables have far too many columns that say roughly the same thing. Using an invented hotel group on the West Cork coast and Northgate Hampers in Cork, this chapter shows how principal component analysis and factor analysis boil them down to a few dimensions you can name and act on, with the maths, working code and the checks that decide whether to trust the result. Part 12 of 21 of the Marketing Analytics series.

Twenty four columns of averages and no answer

Every autumn Carbery Coast Hotels, an invented group of five small hotels in Clonakilty and along the West Cork coast, sends its guests a survey. Twenty four questions, each rated from 1 to 10: how warm the reception was, how quickly staff solved problems, whether the room was clean, how the breakfast tasted, whether the price felt fair. Last season 1,340 guests answered. Áine, who runs the group with her brother, sends me the spreadsheet with a one line email: "Everything is between 7.6 and 8.7. What do we fix first, and will it move direct bookings?"

That is the honest problem with survey data and most customer databases. You do not have 24 pieces of information. You have perhaps four things, each measured six times in slightly different words, plus noise. A guest who gives the reception a 9 gives "staff solved problems quickly" a 9 as well, because both ask how the people treated her. Until you find those underlying things, every chart of 24 bars is the same four bars drawn over and over, and what to fix stays a matter of taste.

Principal component analysis and factor analysis are the two tools for this job. They take dozens of correlated variables and boil them down to a handful of dimensions you can name, score and act on. They look alike, they are often confused, and they answer different questions. This chapter covers both, using Áine's survey and the customer behaviour of Northgate Hampers, Priya's invented Cork hamper shop from earlier chapters. Carbery Coast, Northgate and every number here are made up so I can show each step. The numbers come from simulated files, and the code further down prints exactly what the tables show.

Where this sits in the series

This is chapter 12 of 21 in my Marketing Analytics series, in part three on interrelationship techniques (chapters 11 to 14). The previous chapter, simultaneous equations: when marketing causes sales and sales cause marketing, was the bridge: two outcomes, each also a cause. From here on there is no outcome at all. We stop asking what predicts sales and start asking how variables hang together. The next chapter, segmentation: strategy before algorithms, asks what a segment is for before anyone opens a clustering library, and chapter 14 then uses the components built here as raw material.

Two techniques that look alike and are not

Both methods start from the same object, the correlation matrix of your variables: a square table that says how strongly each pair moves together. Both end with a few new columns that replace the many old ones. The philosophies differ, and the difference matters when you interpret the result.

Principal component analysis (PCA) is geometry. It finds the direction in which the cloud of data points is widest, the first component, then the widest direction at right angles to it, and so on until it has as many new axes as there were variables, describing the same cloud, just turned. Nothing is assumed about why the variables correlate. There is no error term. PCA rewrites the data without losing anything, with the interesting variance pushed to the front so you can drop the tail.

Factor analysis (FA) is a model. It assumes that a few unobserved constructs, the latent factors, cause the observed answers. A guest's invisible opinion of the service, plus some noise specific to each question, produces her ratings of the reception, the staff's speed and their knowledge of the area. Factor analysis splits each variable's variance into a common part explained by the factors and a unique part that belongs to that item alone. PCA keeps all the variance. FA only wants the shared part.1

So when you believe your items measure constructs, as in a survey, factor analysis is the honest tool. When you just want to compress correlated behavioural variables into inputs for clustering or regression, PCA is faster, more stable and assumes less. They often give similar loadings, which is why they are confused, and they diverge when the unique variances are large.

An example: Priya's customer table has an order count and a total spend. Nobody believes a hidden "spend construct" causes both. They move together because spend is mostly orders times basket. That is a PCA job. Áine's three questions about reception, problem solving and local knowledge are different. They were written to measure one thing, the service, from three angles. That is a factor analysis job.

Standardise first, then build the correlation matrix

PCA finds the direction of largest variance in whatever units you hand it. Spend in euro varies in the thousands, a click rate in the hundredths, so on raw data the first component is simply spend. The fix is to put every variable on the same scale first:

zij=xij−xˉjsj,R=1n−1Z⊤Zz_{ij} = \frac{x_{ij} - \bar{x}_j}{s_j}, \qquad R = \frac{1}{n-1} Z^{\top} Z

Here xijx_{ij} is the raw value of variable jj for customer or respondent ii, xˉj\bar{x}_j is the mean of variable jj, sjs_j is its standard deviation, and zijz_{ij} is the standardised value, which says how many standard deviations above or below the average this person sits. ZZ is the table of all standardised values, one row per person and one column per variable, Z⊤Z^{\top} is that table turned on its side (the transpose), nn is the number of rows, and RR is the correlation matrix: a pp by pp table, pp being the number of variables, with ones on the diagonal and the correlation of each pair everywhere else. Means and standard deviations are in Foundations: spread, correlation in Foundations: correlation, and multiplying tables like this in Foundations: matrices.

On raw data the result depends on the units of measurement, so PCA on raw values only makes sense when every variable shares one unit.2

Worked example: two Northgate columns

Take a customer with seven orders when the average Northgate customer has 3.4 and the standard deviation is 2.6. Her standardised order count is (7−3.4)/2.6=1.38(7 - 3.4) / 2.6 = 1.38. Her spend gets the same treatment in its own units. Across 6,200 customers, log orders and log spend correlate at 0.89, so the correlation matrix for just these two variables has ones on the diagonal (each variable with itself) and 0.89 in both other cells. I will carry this small matrix through the next section.

The maths of principal components

The correlation matrix is decomposed into its eigenvectors and eigenvalues:

R=VΛV⊤,Λ=diag(λ1,λ2,…,λp),λ1≥λ2≥⋯≥λp≥0R = V \Lambda V^{\top}, \qquad \Lambda = \mathrm{diag}(\lambda_1, \lambda_2, \dots, \lambda_p), \qquad \lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p \ge 0

Here RR is the pp by pp correlation matrix, VV is a pp by pp matrix whose columns are the eigenvectors (each column holds the weights that define one component), V⊤V^{\top} is its transpose, and Λ\Lambda is a diagonal matrix holding the eigenvalues λ1\lambda_1 to λp\lambda_p, sorted from largest to smallest. Each eigenvalue is the variance of one component.2 The matrix notation is explained in Foundations: matrices, and variance in Foundations: spread.

The kkth component score for a person is a weighted sum of her standardised values:

ck=vk⊤z=v1k z1+v2k z2+⋯+vpk zpc_k = v_k^{\top} z = v_{1k}\, z_1 + v_{2k}\, z_2 + \dots + v_{pk}\, z_p

where ckc_k is her score on component kk, zz is her vector of pp standardised values z1z_1 to zpz_p, vkv_k is the kkth column of VV, and vjkv_{jk} is the weight of variable jj in component kk. A weighted sum of a row of numbers is one row of a matrix product, see Foundations: matrices.

The loading of variable jj on component kk is the correlation between the two, and it is the weight scaled up by the square root of the component's variance:

ℓjk=vjkλk\ell_{jk} = v_{jk} \sqrt{\lambda_k}

Here ℓjk\ell_{jk} is the loading, vjkv_{jk} the weight from above and λk\lambda_k the eigenvalue of component kk. A loading reads like any correlation, from minus 1 to 1, see Foundations: correlation.

Because every standardised variable has a variance of exactly one, the variances on the diagonal of RR add up to pp, and so do the eigenvalues.2 The share of total variance carried by component kk is then simply

sharek=λk∑j=1pλj=λkp\text{share}_k = \frac{\lambda_k}{\sum_{j=1}^{p} \lambda_j} = \frac{\lambda_k}{p}

where sharek\text{share}_k is the proportion of all variance that component kk explains, λk\lambda_k is its eigenvalue, ∑j=1pλj\sum_{j=1}^{p} \lambda_j is the sum of all pp eigenvalues and pp the number of variables. It is an ordinary share of a total, see Foundations: ratios and rates.

Worked example: orders and spend become one component

For a two by two correlation matrix with correlation rr, the eigenvalues are 1+r1 + r and 1−r1 - r, and the eigenvectors are the two diagonals, with weights of 1/2=0.70711/\sqrt{2} = 0.7071 each.

  • Eigenvalues: 1+0.89=1.891 + 0.89 = 1.89 and 1−0.89=0.111 - 0.89 = 0.11. They add up to 2, the number of variables.
  • Share of the first component: 1.89/2=0.9451.89 / 2 = 0.945, so 94.5% of everything orders and spend say, they say together.
  • First component weights: v11=v21=0.7071v_{11} = v_{21} = 0.7071. The component is "orders and spend, equally".
  • Loadings: ℓ11=0.7071×1.89=0.7071×1.375=0.97\ell_{11} = 0.7071 \times \sqrt{1.89} = 0.7071 \times 1.375 = 0.97. Both variables correlate 0.97 with the component.
  • A customer with standardised orders 1.2 and spend 0.8 scores 0.7071×1.2+0.7071×0.8=1.410.7071 \times 1.2 + 0.7071 \times 0.8 = 1.41 on the first component and 0.7071×(1.2−0.8)=0.280.7071 \times (1.2 - 0.8) = 0.28 on the second, the small leftover where she spends a little less per order than her order count suggests.

Two columns, one real dimension. With 18 columns the computer does the arithmetic, and nothing else changes.

How many components or factors to keep

Kaiser's rule keeps every component with an eigenvalue above one, since a component worth less than a single original variable compresses nothing. It is the default in most software, and in Costello and Osborne's simulations it kept too many factors in 36% of samples.1 The scree plot draws the eigenvalues in order and looks for the elbow where the curve flattens into rubble.2 Parallel analysis, proposed by John Horn in 1965, computes the eigenvalues of random data of the same size and keeps only components whose eigenvalues beat what randomness alone produces.3 Costello and Osborne call it accurate and easy to use.

I use all three. When they agree, I trust the count. When they disagree, I keep the smaller number and say so. And interpretability overrules all three: a factor you cannot name from its items is not a finding. Both worked examples below show the three rules at work, and in the Northgate one they disagree.

The factor model

Factor analysis writes each observed variable as a weighted sum of a few latent factors plus its own unique part:

xi=λi1f1+λi2f2+⋯+λimfm+ei,or in matrix formx=Lf+ex_i = \lambda_{i1} f_1 + \lambda_{i2} f_2 + \dots + \lambda_{im} f_m + e_i, \qquad \text{or in matrix form} \qquad \mathbf{x} = L \mathbf{f} + \mathbf{e}

Here xix_i is the standardised score on item ii, f1f_1 to fmf_m are the mm latent factors (far fewer than the pp items), λik\lambda_{ik} is the loading of item ii on factor kk, and eie_i is the unique part of item ii, the bit no factor explains. In matrix form x\mathbf{x} holds all pp items, LL is the pp by mm loading matrix and e\mathbf{e} the unique parts. Each line is a regression on variables nobody can see, see Foundations: regression and Foundations: matrices.

If the factors are uncorrelated with unit variance, the model implies a structure for the whole correlation matrix:4

R=LL⊤+Ψ,hi2=∑k=1mλik2,ψi=1−hi2R = L L^{\top} + \Psi, \qquad h_i^2 = \sum_{k=1}^{m} \lambda_{ik}^2, \qquad \psi_i = 1 - h_i^2

where RR is the correlation matrix, LL⊤L L^{\top} is the part the factors explain, Ψ\Psi is a diagonal matrix of the unique variances ψi\psi_i, and hi2h_i^2 is the communality of item ii: the sum of its squared loadings λik2\lambda_{ik}^2 over the mm factors, which is the share of its variance that the factors explain. Shares of variance are in Foundations: spread and squared correlations in Foundations: correlation.

An item with a communality of 0.70 is well represented. One with 0.25 is mostly noise or measures something the factors do not cover, and should be questioned or dropped. Estimation is by maximum likelihood or principal axis factoring. Both iterate until the implied and observed correlations agree as closely as possible.4

Worked example: the communality of one survey item

In Áine's final solution, "reception was warm and welcoming" loads 0.73 on the first factor, 0.16 on the second, 0.15 on the third and 0.18 on the fourth. Its communality is

h2=0.732+0.162+0.152+0.182=0.533+0.026+0.023+0.032=0.61h^2 = 0.73^2 + 0.16^2 + 0.15^2 + 0.18^2 = 0.533 + 0.026 + 0.023 + 0.032 = 0.61

so the four factors explain 61% of how guests rate the reception, and the unique part is ψ=1−0.61=0.39\psi = 1 - 0.61 = 0.39. That 39% belongs to the item: the particular receptionist on the particular evening, and noise.

Rotation changes the reading, not the fit

A factor solution is not unique: any rotation of the factor axes fits the data exactly as well.4 That is a freedom, not a flaw, and you use it to pick the rotation that is easiest to read. Varimax, introduced by Henry Kaiser in 1958, keeps the factors at right angles and pushes each item to load strongly on one factor and weakly on the rest.5 Oblimin lets the factors correlate, which is more realistic for surveys, where a guest who feels well served also tends to feel the price was fair.

I start with varimax because the table is cleaner, then check oblimin to see whether the story changes. Costello and Osborne go further and recommend oblique rotation by default, because it reproduces an orthogonal solution when the factors really are uncorrelated and keeps the information when they are not.1 Either way, choose the rotation before you see the loadings, and report it.

Worked example: turning the axes by 45 degrees

An item loads 0.60 on factor one and 0.60 on factor two. Its communality is 0.602+0.602=0.720.60^2 + 0.60^2 = 0.72. Turn both axes by 45 degrees, so that cos⁡45∘=sin⁡45∘=0.7071\cos 45^\circ = \sin 45^\circ = 0.7071:

  • New loading on factor one: 0.60×0.7071+0.60×0.7071=0.850.60 \times 0.7071 + 0.60 \times 0.7071 = 0.85.
  • New loading on factor two: 0.60×0.7071−0.60×0.7071=00.60 \times 0.7071 - 0.60 \times 0.7071 = 0.
  • New communality: 0.852+02=0.720.85^2 + 0^2 = 0.72.

The item now clearly belongs to one factor, and its communality has not moved. Rotation buys readability at no cost in fit.

Worked example: Carbery Coast Hotels, 24 items, four factors

Áine's 1,340 responses, 24 items, 24 eigenvalues adding up to 24. The Kaiser Meyer Olkin measure, which asks whether the items share enough variance to factor at all, comes out at 0.94. Kaiser called values above 0.9 "marvelous" and anything below 0.5 unacceptable.6 So the data is fit to factor. The first 12 eigenvalues, against the 95th percentile of the same eigenvalues from random data of the same shape:

Scree plot of the 24 survey items against random dataThe eigenvalues of the 24 Carbery Coast survey items fall from 7.69 to 2.70, 2.04 and 1.85, then drop to 0.71 at the fifth and flatten towards 0.51 at the twelfth. Random data of the same shape gives 1.28 falling slowly to 1.02. The first four real eigenvalues lie above the random line, the fifth and all later ones below it.Eigenvalues of the 24 survey items02.557.510135791112Survey dataRandom data, 95th percentileEigenvalueComponent number

The horizontal axis is the component number, the vertical axis its eigenvalue. The survey line drops steeply to the fourth component and then flattens, while the random line runs almost level just above one. Where the survey line crosses below it, between the fourth and fifth component, structure ends and noise begins.

The first component alone carries 7.69 of 24, which is 32.0% of all variance: a general "everything was fine" tendency running through every answer. The second carries 11.3%, the third 8.5%, the fourth 7.7%. Together the four explain 14.28/24=59.5%14.28 / 24 = 59.5\%. The fifth eigenvalue, 0.71, is below Kaiser's line of one and below the 1.15 random data gives in fifth place, and the scree flattens. All three rules agree: four.

Then the maximum likelihood fit with varimax rotation. One item did not survive it. "WiFi speed" loaded 0.38 on the room factor, 0.36 on the service factor and 0.15 and 0.10 on the other two, for a communality of 0.382+0.362+0.152+0.102=0.310.38^2 + 0.36^2 + 0.15^2 + 0.10^2 = 0.31. It measures none of the four constructs. I take it out of the solution and report it separately, the right treatment for an item that matters operationally but belongs to no construct. The refit on the remaining 23 items gives these loadings, 12 of them shown (the other 11 follow the same pattern):

Survey itemServiceRoomFoodValue
Reception was warm and welcoming0.730.160.150.18
Staff solved problems quickly0.750.170.120.18
Staff knew the area well0.630.150.130.16
Room was clean on arrival0.170.790.110.10
Bed and sleep quality0.180.730.070.07
Bathroom condition0.160.740.110.09
Breakfast quality0.170.110.770.19
Dinner quality0.140.080.720.19
Local West Cork produce on the menu0.160.100.620.11
Price felt fair for what we received0.190.150.150.73
Extras were fairly priced0.160.130.160.69
Would book direct again at this price0.270.140.170.57

Read it row by row. The first three items load between 0.63 and 0.75 on the first factor and below 0.19 on the other three, so factor one is about people, and I call it Service. Rows four to six cluster on the second factor with loadings from 0.73 to 0.79 and nothing above 0.18 elsewhere: Room. Rows seven to nine define Food, with "local West Cork produce on the menu" weaker at 0.62: partly about food, partly about whether the guest cares about local produce at all. The last three rows form Value. Note the last one: "would book direct again at this price" loads 0.57 on Value but also 0.27 on Service, more than any other item's second loading. That mild cross loading is informative: the decision to come back is mostly about fairness and a little about how you were treated.

The structure the model assumes, latent causes first and observed answers after:

The factor model behind the Carbery Coast surveyFour unseen opinions, Service, Room, Food and Value, each drive their own group of ratings. Every rating also carries a unique part of its own. Together they produce the 23 ratings in the spreadsheet, and the WiFi rating sits outside all four.Four unseen opinions: Service, Room,Food, ValueService drives reception, problemsolving, local knowledge and three moreratingsRoom drives cleanliness, sleep,bathroom and three moreFood drives breakfast, dinner, localproduce and two moreValue drives fair price, fair extras,booking direct again and three moreEach rating adds its own unique part:the evening, the receptionist, noiseWiFi speed belongs to none of the fourand is reported on its ownThe 23 ratings you see in thespreadsheet

The arrows run from the invisible to the visible. The analysis works the other way, from the 23 ratings back to the four opinions that most plausibly produced them.

From four factors to one decision

Naming factors is pleasant, but it is not a decision. The survey also asked "how likely are you to recommend us" on a 0 to 10 scale. So I compute each guest's four factor scores (centred on zero, one unit roughly one standard deviation) and regress the recommendation score on them:

PredictorCoefficientStandard errorReading
Intercept8.140.03the average guest at average factor scores
Service0.560.03one unit more service adds 0.56 points
Room0.440.03the smallest lever of the four
Food0.530.03close behind Service
Value0.690.03the largest lever, and it is not the room rate

The model explains 48% of the variance in recommendation. Because the varimax factors are uncorrelated, their estimated scores are nearly so (no pair above 0.07), and the coefficients do not fight each other the way 24 raw items would. You can read each one on its own. Regression and its standard errors are covered in modelling demand and price elasticity with regression.

Value is the strongest lever, 0.69/0.44=1.60.69 / 0.44 = 1.6 times the effect of Room. But the Value factor is built as much from "extras were fairly priced" as from the room rate. In the example's follow up survey, which splits the extras question into its parts, parking, spa access and the minibar are the sore points: pricing policy decisions that cost almost nothing to change, not a discount on the room. Room, where Áine had planned to spend most of next year's budget, is the weakest of the four levers for the guests she already has.

That is the whole payoff: 24 columns collapsed into four, ranked against the one number the business cares about. Whether guests who book through the big platforms react differently is a question for the next chapter, which is why the factor scores stay as columns in the guest table instead of dying with the report.

Worked example: Northgate Hampers, 18 variables, three components

The second example needs no survey. For each of 6,200 Northgate customers with an order in the last 14 months, Priya's order history gives 18 behavioural variables: how much and how often they buy, how they gift (message cards, other delivery addresses, December, a company on the invoice), how they respond to offers (vouchers, discount depth, email clicks, days from offer to order), and four that turn out to matter less (orders placed on a phone, newsletter, reviews, returns).

Chapter 14 wants to cluster these customers, and feeding it 18 correlated columns would be a mistake. Orders and spend correlate at 0.89, so k means would count the same information twice, and the "segments" would be spend bands in disguise. This is a PCA job: no theory that a latent "gifting orientation" causes anything, just the fewest uncorrelated columns that keep most of the information. Five skewed variables are logged first, because a few corporate buyers spend many times the median and would otherwise make the first component "the whales against everybody else". Then standardise and decompose:

Scree plot of the 18 Northgate variables against random dataThe eigenvalues of the 18 Northgate behavioural variables are 3.55, 2.65 and 2.59, then drop to 1.05, 1.02, 1.00, 1.00, 0.96, 0.77 and 0.69. Random data of the same shape gives 1.10 falling slowly to 1.00. The first three real eigenvalues lie far above the random line, the fourth sits just below it.Eigenvalues of the 18 behavioural variables01.32.53.851357910Northgate dataRandom data, 95th percentileEigenvalueComponent number

The horizontal axis is the component number, the vertical axis the eigenvalue. Three components stand well clear of the random line, and from the fourth on the Northgate line runs along it, at or just below. The fourth eigenvalue, 1.05, passes Kaiser's rule and fails parallel analysis, which gives 1.06 in fourth place. Components four to eight, all close to 1.0, are leftovers: what average order value does not share with the rest, and the four variables that correlate with almost nothing (phone share, newsletter, reviews, returns), shuffled among themselves. I keep three, which carry (3.55+2.65+2.59)/18=8.79/18=48.9%(3.55 + 2.65 + 2.59) / 18 = 8.79 / 18 = 48.9\% of the variance. After a varimax rotation the components are easy to name:

ComponentShare of variance after rotationVariables with the largest loadings
1 Volume19.4%orders 0.95, spend 0.90, distinct categories 0.85, months active 0.81, days since last order minus 0.61
2 Gifting orientation14.9%gift message share 0.87, ship to other address 0.86, December share 0.81, company on invoice 0.55, average order value 0.43
3 Deal sensitivity14.5%voucher share 0.87, average discount depth 0.85, offer click rate 0.79, days from offer to order minus 0.70

Rotation moves variance between the three components but not out of them: 19.4+14.9+14.5=48.8%19.4 + 14.9 + 14.5 = 48.8\%, the same total as before rotation give or take rounding.

The sign matters here. "Days since last order" loads minus 0.61 on Volume because a customer who buys a lot was, almost by definition, here recently. A negative loading says the variable moves against the component. "Days from offer email to order" loads negatively on Deal sensitivity for the same reason: the deal sensitive customer orders the day the voucher lands.

Component two is where the shop's real structure sits. It separates the gift buyer, sending hampers to other addresses with a message card and spending a little more per order, from the self buyer with the same order count who ships everything home. Nobody told it what a gift buyer is, a good sign that the distinction lives in the data and not just in Priya's head.

Three uncorrelated scores per customer now go into chapter 14, so k means treats each dimension fairly instead of triple counting spend. Why behaviour beats demographics as a basis for segments is argued in segmentation is not demographics. This is the mechanical step that makes it possible.

Running it yourself

You need a wide table, one row per respondent or customer, one column per variable, with few gaps. For survey data, a common rule of thumb asks for ten responses per item. In Costello and Osborne's simulations even that ratio recovered the correct factor structure only 60% of the time, against 70% at 20 responses per item and 40% at five.1 So 24 items with 1,340 answers, about 56 per item, is comfortable, and 24 items with 90 answers is a coin toss at best. Behavioural variables are usually heavily skewed, so log counts and amounts before standardising.

In Python the PCA side for Northgate is short. It uses scikit learn for the decomposition7 and the factor_analyzer package for the varimax rotation:8

import numpy as np, pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from factor_analyzer.rotator import Rotator

X = pd.read_csv("hamper_behaviour.csv").set_index("customer_id")
skewed = ["orders", "spend", "avg_order_value", "days_since_last", "days_offer_to_order"]
X[skewed] = np.log1p(X[skewed])                  # tame the long tails first
Z = StandardScaler().fit_transform(X)            # PCA on the correlation matrix

eig = np.linalg.eigvalsh(np.corrcoef(Z, rowvar=False))[::-1]       # largest first
rng = np.random.default_rng(1)                   # parallel analysis: same shape, pure noise
rand = np.percentile([np.linalg.eigvalsh(np.corrcoef(rng.standard_normal(Z.shape), rowvar=False))[::-1]
                      for _ in range(50)], 95, axis=0)
print(np.round(eig[:6], 2), np.round(rand[:6], 2), np.round(eig[:3].sum() / len(eig), 3))

pca = PCA(n_components=3).fit(Z)
loadings = pca.components_.T * np.sqrt(eig[:3])  # correlation of each variable with each component
rotated = pd.DataFrame(Rotator(method="varimax").fit_transform(loadings),
                       index=X.columns, columns=["volume", "gifting", "deal"])   # named after reading them
print(rotated.round(2))
scores = pd.DataFrame(Z @ np.linalg.pinv(rotated.values.T), index=X.index, columns=rotated.columns)
scores.to_csv("hamper_components.csv")           # three columns for the segmentation chapter

It prints the eigenvalues and random thresholds from the chart, the share kept (0.489) and the loadings in the table. The three saved score columns correlate with each other at exactly zero.

For the survey, the same factor_analyzer package fits the maximum likelihood factor model with rotation and computes the Kaiser Meyer Olkin measure:8

import numpy as np, pandas as pd
import statsmodels.api as sm
from factor_analyzer import FactorAnalyzer
from factor_analyzer.factor_analyzer import calculate_kmo

survey = pd.read_csv("guest_survey.csv")
items = survey.drop(columns=["guest_id", "recommend"])        # the 24 rating items
print("KMO", round(calculate_kmo(items)[1], 2))

eig = np.linalg.eigvalsh(items.corr())[::-1]
rng = np.random.default_rng(1)                   # parallel analysis, as for the hamper shop
rand = np.percentile([np.linalg.eigvalsh(np.corrcoef(rng.standard_normal(items.shape), rowvar=False))[::-1]
                      for _ in range(50)], 95, axis=0)
print(np.round(eig[:6], 2), np.round(rand[:6], 2))

fit = FactorAnalyzer(n_factors=4, rotation="varimax", method="ml").fit(items)
h2 = pd.Series(fit.get_communalities(), index=items.columns)
print(h2[h2 < 0.4].round(2))                     # items the four factors barely explain
items = items.drop(columns=h2[h2 < 0.4].index)
fit = FactorAnalyzer(n_factors=4, rotation="varimax", method="ml").fit(items)
print(pd.DataFrame(fit.loadings_, index=items.columns).round(2))   # the loadings table

names = ["room", "service", "value", "food"]          # named after reading the loadings above
scores = pd.DataFrame(fit.transform(items), columns=names)
model = sm.OLS(survey["recommend"], sm.add_constant(scores)).fit()
print(model.params.round(2), model.bse.round(2), round(model.rsquared, 2))

It prints KMO 0.94, the eigenvalues from the chart, flags wifi_speed 0.31, then the loadings and the regression from the tables, with an R squared of 0.48. In R, the psych package does the same with fa.parallel for the count and fa for the rotated fit.9

Time: a survey factor analysis is a week of work including cleaning, rotation checks and write up, and most of the week is the cleaning. A behavioural PCA is two days once the customer table exists, and the customer table is where the real effort goes.

Five checks before you trust a solution

  1. Sampling adequacy. The Kaiser Meyer Olkin measure should be well above 0.6. Kaiser labelled the 0.60s "mediocre" and the 0.70s "middling".6 Below that, the variables do not share enough variance to factor.
  2. Reliability. Cronbach's alpha for the items of each factor, which says how well the items agree with each other,10 should be at least 0.70, the level usually treated as acceptable in social science research.11
  3. The count. Parallel analysis should confirm the number of factors.
  4. Stability. Split the sample in half at random, fit both halves and compare the loadings. If the halves disagree about which items belong where, the structure is not stable, and I report fewer factors.
  5. Rotation. Check whether an oblique rotation changes the story. If oblimin says two factors correlate at 0.7, they may be one factor, and I should say so rather than present a tidy four when the data prefers three.

An example: Carbery Coast passes the first and third checks above, a KMO of 0.94 and four factors confirmed by parallel analysis. The other three need the item groups, two halves of 670 guests and an oblimin refit.

Pitfalls

The naming fallacy. The moment you call a factor "Value" you begin to believe the data measured value. It did not. It measured six questions that a rotation grouped together, and your name is a hypothesis about why. I once saw a factor labelled "brand love" built entirely from items about delivery speed. Name factors from their items, keep the items next to the name in the report, and treat the name as a working title.

Kaiser's rule on autopilot. It is among the least accurate ways to choose the number of factors and tends to keep too many.1 At Northgate, a handful of loosely related variables produce eigenvalues near one, and Kaiser would happily keep the first of them as a "component".

The agreeable guest. If some guests rate everything highly out of politeness, part of the large first unrotated component is a mood rather than a construct. Rotation helps, reverse worded items help more, and reporting that component as "overall satisfaction" without saying so is a small dishonesty.

Treating components as permanent. A component summarises a particular dataset in a particular season. Rerun the Northgate analysis next year with new products and the gifting component may absorb the corporate variables or split in two. Components are tools for the decision in front of you, not facts about your customers, and dashboards built on them should carry a date.

How I do this for clients

The deliverable depends on which of the two problems you have. If you have a survey, or want one, I need the item wording, the raw responses one row per respondent and, ideally, one business outcome per respondent, such as whether they rebooked or renewed. If you have a customer database, I need the order and event history so that I can build the behavioural table myself, because the variables already sitting in the CRM are rarely the right ones.

It starts with the free workshop, a couple of hours in which we agree the decision the analysis has to support and what the output must look like for you to act on it. Then two weeks of real work on me, before you commit to anything. For a survey that means the factor solution with loadings and reliabilities, factor scores written back to your respondent table, and a two page memo that ranks the factors against your outcome with the coefficient, its error margin and the risk I see in acting on it. For a database it means the behavioural table, the component solution and the scores stored as columns in your own tables, ready for segmentation, plus a memo saying in plain words what each component means and what it does not.

Your raw data stays yours, and the terms for code and models are in the analytics and data science terms. The Data Science page describes the range of work, and the fees are on the pricing page. One accountable person, honest numbers with their uncertainty attached, no vanity metrics.

Questions to ask before you act on a factor solution

  • How many respondents or customers per variable went into this, and were there enough to trust the loadings?
  • Which rule decided the number of factors, and would parallel analysis or the scree plot have given a different count?
  • Were the variables standardised or logged before the decomposition, and which ones dominated if not?
  • Which rotation was used, was it chosen before the loadings were seen, and does an oblique rotation change the story?
  • Do the factor names follow from the items in the table, or from what we hoped to find?
  • Have the factor scores been linked to a business outcome, and how much of that outcome do they explain?
  • Are the scores stored where the next analysis can use them, with a date on them?

The series is inspired by Mike Grigsby's Marketing Analytics.12 The explanations, examples and numbers here are my own.

Footnotes

  1. Anna B. Costello and Jason Osborne, Best practices in exploratory factor analysis: four recommendations for getting the most from your analysis, Practical Assessment, Research, and Evaluation 10(7), 2005, https://openpublishing.library.umass.edu/pare/article/id/1650/

  2. Penn State Eberly College of Science, STAT 505 Applied Multivariate Statistical Analysis, Lesson 11 Principal Components Analysis, https://online.stat.psu.edu/stat505/Lesson11

  3. John L. Horn, A rationale and test for the number of factors in factor analysis, Psychometrika 30, 1965, https://doi.org/10.1007/BF02289447

  4. Penn State Eberly College of Science, STAT 505 Applied Multivariate Statistical Analysis, Lesson 12 Factor Analysis, https://online.stat.psu.edu/stat505/Lesson12

  5. Henry F. Kaiser, The varimax criterion for analytic rotation in factor analysis, Psychometrika 23, 1958, https://doi.org/10.1007/BF02289233

  6. William Revelle, psych package documentation, KMO: Kaiser Meyer Olkin measure of sampling adequacy, CRAN, https://search.r-project.org/CRAN/refmans/psych/html/KMO.html

  7. scikit learn developers, sklearn.decomposition.PCA, scikit learn documentation, https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html

  8. Educational Testing Service, factor_analyzer documentation (FactorAnalyzer, Rotator, calculate_kmo), https://factor-analyzer.readthedocs.io/en/latest/factor_analyzer.html

  9. William Revelle, psych package documentation, fa.parallel: scree plots of data or correlation matrix compared to random data, CRAN, https://search.r-project.org/CRAN/refmans/psych/html/fa.parallel.html

  10. Lee J. Cronbach, Coefficient alpha and the internal structure of tests, Psychometrika 16, 1951, https://doi.org/10.1007/BF02310555

  11. UCLA Office of Advanced Research Computing, What does Cronbach's alpha mean, https://stats.oarc.ucla.edu/spss/faq/what-does-cronbachs-alpha-mean/

  12. Mike Grigsby, Marketing Analytics: A Practical Guide to Improving Consumer Insights Using Data Techniques, Kogan Page, 2022, ISBN 9781398608191, https://www.netgalley.co.uk/catalog/book/267990

If you have a survey with too many questions, or a customer table with too many columns, and nobody can tell you which four things actually matter, bring it to the free workshop. We agree the decision the analysis has to support, then the first two weeks of work are on me, and I tell you honestly whether a factor solution is in there or the data is telling a different story.