import numpy as np, pandas as pd, warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
from scipy.stats import pearsonr
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.model_selection import cross_val_score, KFold
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False,
"axes.grid":True,"grid.alpha":0.22,"axes.titleweight":"bold","axes.titlesize":12.5,"axes.titlelocation":"left"})
VI, DK, LT, MUT, GD, RD = "#6d28d9", "#4c1d95", "#a78bfa", "#94a3b8", "#047857", "#dc2626"
BASE = "https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
fn = "capstone-pca-survey-battery.xlsx"
def load(sheet):
try: return pd.read_excel("../../data/" + fn, sheet_name=sheet)
except FileNotFoundError: return pd.read_excel(BASE + fn, sheet_name=sheet)
raw = load("Data")
battery = load("Battery")
notes = load("Notes")
Q = [c for c in raw.columns if c.startswith("q") and c[1].isdigit()]
print(f"rows {len(raw):,} battery items {len(Q)}")
raw.head()
rows 1,896 battery items 20
| respondent_id | q01_app_easy | q02_site_fast | q03_find_info | q04_checkout_simple | q05_app_crashes | q06_account_setup | q07_app_to_support | q08_agent_knowledge | q09_agent_courtesy | ... | q16_worth_paying | q17_fees_unclear | q18_plan_choice | q19_value_vs_rivals | q20_ads_appealing | region | age_band | tenure_years | overall_satisfaction | would_recommend | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 91042 | 4 | 4 | 4 | 5 | 3 | 5 | 4 | 4 | 4 | ... | 5 | 3 | 4 | 4 | 4 | West | 55+ | 4.3 | 6 | 1 |
| 1 | 90975 | 4 | 4 | 3 | 4 | 3 | 4 | 4 | 5 | 5 | ... | 2 | 5 | 4 | 3 | 1 | Midwest | 35-54 | 2.2 | 6 | 1 |
| 2 | 91408 | 4 | 4 | 4 | 4 | 2 | 5 | 5 | 4 | 5 | ... | 2 | 4 | 2 | 4 | 3 | Northeast | 55+ | 3.1 | 5 | 1 |
| 3 | 90075 | 3 | 2 | 2 | 2 | 4 | 5 | 4 | 4 | 3 | ... | 5 | 2 | 5 | 4 | 5 | South | 55+ | 1.3 | 5 | 1 |
| 4 | 90131 | 3 | 3 | 2 | 3 | 3 | 1 | 5 | 4 | 5 | ... | 4 | 3 | 3 | 2 | 5 | West | 35-54 | 6.1 | 5 | 1 |
5 rows × 26 columns
Step 1 · The brief, and the battery¶
This is the customer survey from Chapter 177, fielded again. This wave carried the twenty-item agreement battery the research team had been arguing about for a year.
for line in notes.Notes.fillna(""):
print(line)
CAPSTONE 32 - TWENTY ITEMS, HOW MANY INDICES? The customer survey from Chapter 177, fielded again. This wave carried the twenty item agreement battery the research team had been arguing about for a year. THE REQUEST. Marketing wants 'a satisfaction score' for the dashboard. The research team thinks there are three separate things in here and that averaging them into one number would hide the one that is moving. Nobody has looked at the correlations. THE JOB. Decide how many distinct things the battery measures, decide which items belong to which, name them honestly, and produce indices the dashboard can compute every month without re-running an analysis. SCALE. Every battery item is agreement, 1 = strongly disagree to 5 = strongly agree. FOUR ITEMS ARE NEGATIVELY WORDED (q05, q11, q13, q17). They are stored exactly as collected. Agreeing with 'I wait too long to reach someone' is a BAD outcome, and agreeing with 'staff treat me courteously' is a good one, on the same 1 to 5 scale. KNOWN FAULTS IN THIS EXPORT, left in deliberately: - a block of 46 responses was exported twice - 58 respondents gave the identical answer to all twenty items - 'prefer not to say' was coded as 99 rather than left blank WHAT THE DASHBOARD ACTUALLY NEEDS. Whatever comes out of this has to be recomputable in the reporting tool every month, by someone who will not have this notebook.
pd.set_option("display.max_colwidth", 70)
print(battery.to_string(index=False))
column statement
q01_app_easy The app is easy to use.
q02_site_fast The website loads quickly.
q03_find_info I can find the information I need without help.
q04_checkout_simple Paying and checking out is simple.
q05_app_crashes The app crashes or freezes on me. [NEGATIVELY WORDED]
q06_account_setup Setting up my account was straightforward.
q07_app_to_support It is easy to reach a person from inside the app.
q08_agent_knowledge The agent I spoke to knew what they were doing.
q09_agent_courtesy Staff treat me courteously.
q10_first_contact My issue was handled the first time I raised it.
q11_wait_too_long I wait too long to reach someone. [NEGATIVELY WORDED]
q12_kept_informed I am kept informed while an issue is open.
q13_repeat_myself I have to repeat myself to different people. [NEGATIVELY WORDED]
q14_issue_resolved My issues get properly resolved.
q15_price_fair The price I pay is fair.
q16_worth_paying The service is worth what it costs.
q17_fees_unclear The fees are hard to understand. [NEGATIVELY WORDED]
q18_plan_choice There is a plan that suits how I use the service.
q19_value_vs_rivals Compared with alternatives, this is good value.
q20_ads_appealing I enjoy the company's advertising.
Note the four items marked negatively worded. They sit on the same 1-to-5 agreement scale as the others, which means agreeing with "I wait too long to reach someone" is a bad outcome while agreeing with "staff treat me courteously" is a good one. Section 4 is about what happens if that is forgotten, and the answer is not what most people expect.
Step 2 · Cleaning¶
d = raw.drop_duplicates(subset="respondent_id").copy()
print(f"raw rows {len(raw):,}")
print(f"after the duplicates {len(d):,}")
d[Q] = d[Q].replace(99, np.nan) # "prefer not to say" was coded 99, not left blank
d = d.dropna(subset=Q)
print(f"complete on all 20 {len(d):,}")
flat = d[Q].std(axis=1) == 0 # identical answer to every item
print(f"straight-liners {int(flat.sum())} (removed)")
d = d[~flat].reset_index(drop=True)
print(f"analysis sample {len(d):,}")
raw rows 1,896 after the duplicates 1,850 complete on all 20 1,677 straight-liners 49 (removed) analysis sample 1,628
Straight-liners matter more here than in most analyses. A respondent who answered 4 to all twenty items contributes a row of zero variance that is perfectly consistent with every other straight-liner, and correlations are exactly what this whole chapter is built on.
The 99 codes matter for the same reason. Left in, "prefer not to say" would enter every correlation as a value nineteen points above the top of the scale.
Step 3 · First look¶
Component analysis is an analysis of a correlation matrix. It is worth looking at the matrix before handing it to an algorithm, because everything the method can possibly find is already in this picture.
R = d[Q].astype(float).corr()
off = R.values[~np.eye(len(Q), dtype=bool)]
avg = ((R.sum() - 1)/(len(Q) - 1)).sort_values()
neg = avg[avg < 0]
print(f"{len(Q)} items, {len(d):,} respondents, {len(off)//2} distinct correlations")
print(f" strongest {off.max():+.3f} positive and {off.min():+.3f} negative; "
f"{(off < 0).mean():.1%} of all pairs are negative")
print(f" items whose average correlation with the other {len(Q)-1} is NEGATIVE "
f"({len(neg)} of {len(Q)}):")
print(" " + ", ".join(f"{c} ({v:+.2f})" for c, v in neg.items()))
print(f" responses average {d[Q].values.mean():.2f} on the 1 to 5 scale "
f"(sd {d[Q].values.std():.2f}); lowest {d[Q].mean().idxmin()} at {d[Q].mean().min():.2f}, "
f"highest {d[Q].mean().idxmax()} at {d[Q].mean().max():.2f}")
20 items, 1,628 respondents, 190 distinct correlations
strongest +0.579 positive and -0.525 negative; 33.7% of all pairs are negative
items whose average correlation with the other 19 is NEGATIVE (4 of 20):
q05_app_crashes (-0.20), q11_wait_too_long (-0.20), q13_repeat_myself (-0.17), q17_fees_unclear (-0.17)
responses average 3.33 on the 1 to 5 scale (sd 1.15); lowest q05_app_crashes at 2.81, highest q09_agent_courtesy at 3.65
What this shows. Four items correlate negatively with almost everything else in the battery. Nothing is wrong with those four. They are the ones worded so that agreeing is bad news, and until that is undone every method below will read them as a separate dimension of the survey rather than as part of the thing they measure. Finding them by looking rather than by being told is the point, because on a real battery nobody hands you the list.
One other item is worth a note. Ads appealing sits at +0.05, near zero against everything. It is not reversed, it is unrelated, and it is measuring something the rest of this battery is not.
R = d[Q].astype(float).corr()
short = [c.split("_", 1)[1].replace("_", " ") for c in Q]
fig, axes = plt.subplots(1, 2, figsize=(12.6, 5.0))
ax = axes[0]
im = ax.imshow(R.values, cmap="RdBu_r", vmin=-0.7, vmax=0.7)
ax.set_xticks(range(len(Q))); ax.set_xticklabels([c[:3] for c in Q], fontsize=7.5, rotation=90)
ax.set_yticks(range(len(Q))); ax.set_yticklabels([c[:3] for c in Q], fontsize=7.5)
ax.set_title("Twenty items, and four that run the other way")
ax.grid(False)
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.03).ax.tick_params(labelsize=8)
ax = axes[1]
avg = ((R.sum() - 1) / (len(Q) - 1))
order = np.argsort(avg.values)
cols = [RD if avg.values[i] < 0 else VI for i in order]
ax.barh(range(len(Q)), avg.values[order], color=cols, height=0.68)
ax.set_yticks(range(len(Q)))
ax.set_yticklabels([short[i] for i in order], fontsize=8.5)
ax.axvline(0, color=DK, lw=1.6)
ax.set_xlabel("average correlation with the other nineteen items")
ax.set_title("The four that need attention before anything else")
ax.grid(axis="y", alpha=0)
plt.tight_layout(); plt.show()
Left: three warm blocks sit along the diagonal, which is the structure the battery was designed to have, crossed by four cool stripes and one nearly white row at the bottom. Right: the same four items, isolated, and Step 4 undoes them. Notice what the picture does not say. It shows that items group together, and it says nothing about how many groups there are. That question takes the next three steps and it does not have a single right answer.
Step 4 · Reverse-score the four negatively worded items¶
REV = ["q05_app_crashes", "q11_wait_too_long", "q13_repeat_myself", "q17_fees_unclear"]
e = d.copy()
for c in REV:
e[c] = 6 - e[c] # 1 becomes 5, 5 becomes 1, on a 1 to 5 scale
print("before and after, first five respondents:")
print(pd.concat([d[REV].head(5).add_suffix(" (raw)"),
e[REV].head(5).add_suffix(" (rev)")], axis=1).to_string(index=False))
before and after, first five respondents:
q05_app_crashes (raw) q11_wait_too_long (raw) q13_repeat_myself (raw) q17_fees_unclear (raw) q05_app_crashes (rev) q11_wait_too_long (rev) q13_repeat_myself (rev) q17_fees_unclear (rev)
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 2.0 2.0 5.0 3.0 4.0 4.0 1.0
2.0 2.0 3.0 4.0 4.0 4.0 3.0 2.0
4.0 3.0 4.0 2.0 2.0 3.0 2.0 4.0
3.0 2.0 4.0 3.0 3.0 4.0 2.0 3.0
Step 5 · How many things is this battery measuring?¶
Three rules are in common use. They do not agree, and one of them is much worse than its popularity suggests.
def loadings(X):
Z = (X - X.mean())/X.std(ddof=0)
p = PCA().fit(Z)
return p, p.components_.T*np.sqrt(p.explained_variance_), Z
p1, L1, Z = loadings(e[Q].astype(float))
ev = p1.explained_variance_
print("eigenvalues:", np.round(ev[:8], 3))
print("share of the twenty items' variance:", np.round(p1.explained_variance_ratio_[:8]*100, 1))
print()
print(f"KAISER, keep every eigenvalue above 1: keeps {(ev > 1).sum()}")
print(f" the fourth eigenvalue is {ev[3]:.3f}, which is {1-ev[3]:.3f} below the cutoff")
print()
print("SCREE, look for the elbow: 6.17, 2.56, 2.18, 0.95, 0.66 ...")
print(" the drop from 1 to 2 is large and the drop from 3 to 4 is large.")
print(" Read it as two, or read it as four. This is why the rule is criticized.")
eigenvalues: [6.174 2.559 2.178 0.952 0.659 0.619 0.6 0.577] share of the twenty items' variance: [30.9 12.8 10.9 4.8 3.3 3.1 3. 2.9] KAISER, keep every eigenvalue above 1: keeps 3 the fourth eigenvalue is 0.952, which is 0.048 below the cutoff SCREE, look for the elbow: 6.17, 2.56, 2.18, 0.95, 0.66 ... the drop from 1 to 2 is large and the drop from 3 to 4 is large. Read it as two, or read it as four. This is why the rule is criticized.
Now the third rule, and the one that should be the default. Parallel analysis asks how large an eigenvalue this dataset would produce if there were no structure at all. Simulate datasets of the same shape from uncorrelated variables, take the distribution of each eigenvalue, and keep only the components that beat what noise produces.
rs = np.random.default_rng(4)
sim = np.array([np.linalg.eigvalsh(np.corrcoef(rs.normal(size=(len(e), len(Q))), rowvar=False))[::-1]
for _ in range(500)])
p95 = np.percentile(sim, 95, axis=0)
print("500 random datasets of the same shape")
print(f" {'component':>10s} {'actual':>9s} {'noise 95th':>12s} {'keep?':>7s}")
for i in range(6):
print(f" {i+1:10d} {ev[i]:9.3f} {p95[i]:12.3f} {'yes' if ev[i] > p95[i] else 'no':>7s}")
K = int((ev > p95).sum())
print(f"\nparallel analysis keeps {K}")
print(f" it rejects the fourth component by {p95[3]-ev[3]:.3f}, where Kaiser was within {1-ev[3]:.3f} of keeping it")
500 random datasets of the same shape
component actual noise 95th keep?
1 6.174 1.235 yes
2 2.559 1.189 yes
3 2.178 1.159 yes
4 0.952 1.135 no
5 0.659 1.112 no
6 0.619 1.092 no
parallel analysis keeps 3
it rejects the fourth component by 0.183, where Kaiser was within 0.048 of keeping it
Here all three rules can be made to give the same answer, so the disagreement is not obvious. Apply Kaiser's rule to data with no structure whatever and it becomes obvious.
rs2 = np.random.default_rng(9)
cnt = [(np.linalg.eigvalsh(np.corrcoef(rs2.normal(size=(len(e), 20)), rowvar=False)) > 1).sum()
for _ in range(400)]
print(f"twenty completely uncorrelated variables, n = {len(e):,}, 400 replications")
print(f" 'eigenvalue greater than 1' retains {np.mean(cnt):.1f} components on average")
print(f" range across replications: {min(cnt)} to {max(cnt)}")
print("\nThere is nothing there. The rule finds ten components in it.")
twenty completely uncorrelated variables, n = 1,628, 400 replications 'eigenvalue greater than 1' retains 9.8 components on average range across replications: 8 to 12 There is nothing there. The rule finds ten components in it.
The reason is easy to state. Eigenvalues of a sample correlation matrix scatter around 1 even when the population correlation matrix is the identity, so about half of them land above it. A threshold of exactly 1 has no sampling theory behind it; it is the average of the eigenvalues, not a test.
Parallel analysis is the same idea done properly: compare against what noise actually produces at this sample size and this number of items, rather than against a constant.
Three components. They account for 54.6 percent of the variance in the twenty items.
Step 6 · Rotation, and what it does not do¶
print("UNROTATED LOADINGS, first three components")
print(f" {'item':22s}" + "".join(f"{'PC'+str(j+1):>9s}" for j in range(3)))
for i, c in enumerate(Q):
print(f" {c:22s}" + "".join(f"{L1[i,j]:9.2f}" for j in range(3)))
print(f"\nitems loading above 0.30 on PC1: {(np.abs(L1[:,0]) > 0.30).sum()} of 20")
UNROTATED LOADINGS, first three components item PC1 PC2 PC3 q01_app_easy 0.59 -0.38 -0.36 q02_site_fast 0.58 -0.37 -0.35 q03_find_info 0.54 -0.34 -0.35 q04_checkout_simple 0.62 -0.35 -0.32 q05_app_crashes 0.53 -0.34 -0.29 q06_account_setup 0.49 -0.37 -0.30 q07_app_to_support 0.74 0.05 -0.22 q08_agent_knowledge 0.64 0.49 0.02 q09_agent_courtesy 0.63 0.46 -0.00 q10_first_contact 0.59 0.43 -0.01 q11_wait_too_long 0.58 0.45 0.00 q12_kept_informed 0.55 0.45 0.02 q13_repeat_myself 0.54 0.40 0.01 q14_issue_resolved 0.61 0.46 -0.04 q15_price_fair 0.49 -0.24 0.62 q16_worth_paying 0.47 -0.27 0.60 q17_fees_unclear 0.44 -0.28 0.55 q18_plan_choice 0.42 -0.26 0.56 q19_value_vs_rivals 0.66 -0.32 0.26 q20_ads_appealing 0.17 -0.05 0.23 items loading above 0.30 on PC1: 19 of 20
That is not a description of three constructs. Nineteen of the twenty items load on the first component, because the first component of a correlation matrix with all-positive correlations is always a general "how happy is this person overall" factor. It is mathematically optimal, in the sense that no other three directions capture more variance, and it is not readable.
Rotation fixes the readability by choosing a different basis for the same three-dimensional space.
def varimax(L, tol=1e-7, it=200):
p, k = L.shape; R = np.eye(k); dold = 0
for _ in range(it):
Lam = L @ R
u, s, vt = np.linalg.svd(L.T @ (Lam**3 - Lam @ np.diag((Lam**2).sum(0))/p))
R = u @ vt; dnew = s.sum()
if dnew < dold*(1 + tol): break
dold = dnew
return L @ R
LV = varimax(L1[:, :K].copy())
LV = LV * np.sign(LV.sum(0))
print(f"variance explained, unrotated: {ev[:K].sum()/len(Q)*100:.1f}%")
print(f"variance explained, rotated : {(LV**2).sum()/len(Q)*100:.1f}%")
print("\nThe same. Rotation cannot change how much variance three directions capture.")
variance explained, unrotated: 54.6% variance explained, rotated : 54.6% The same. Rotation cannot change how much variance three directions capture.
Rotation is a choice of coordinates, not a discovery. The three-dimensional subspace is fixed by the data; varimax picks the axes inside it that put each item's weight mostly on one axis. Different criteria pick different axes, all equally valid as summaries and differently readable.
This is worth being clear about with an audience, because "we rotated the factors" often gets heard as "we found the real ones". Nothing was found. The picture was turned around.
Step 7 · The loadings, and reading them honestly¶
print("ROTATED LOADINGS (blank below 0.30)")
print(f" {'item':22s}" + "".join(f"{'C'+str(j+1):>9s}" for j in range(K)) + f"{'communality':>13s}")
for i, c in enumerate(Q):
row = "".join(f"{LV[i,j]:9.2f}" if abs(LV[i,j]) >= 0.30 else f"{'':>9s}" for j in range(K))
print(f" {c:22s}{row}{(LV[i]**2).sum():13.2f}")
multi = (np.abs(LV) > 0.30).sum(1) > 1
print(f"\nitems loading on more than one component: {multi.sum()} of 20 ({', '.join(np.array(Q)[multi])})")
ROTATED LOADINGS (blank below 0.30) item C1 C2 C3 communality q01_app_easy 0.77 0.62 q02_site_fast 0.76 0.60 q03_find_info 0.71 0.53 q04_checkout_simple 0.75 0.60 q05_app_crashes 0.68 0.49 q06_account_setup 0.68 0.47 q07_app_to_support 0.55 0.53 0.60 q08_agent_knowledge 0.79 0.64 q09_agent_courtesy 0.76 0.61 q10_first_contact 0.71 0.54 q11_wait_too_long 0.72 0.54 q12_kept_informed 0.71 0.51 q13_repeat_myself 0.65 0.45 q14_issue_resolved 0.75 0.58 q15_price_fair 0.81 0.68 q16_worth_paying 0.80 0.65 q17_fees_unclear 0.74 0.57 q18_plan_choice 0.74 0.56 q19_value_vs_rivals 0.43 0.62 0.60 q20_ads_appealing 0.08 items loading on more than one component: 2 of 20 (q07_app_to_support, q19_value_vs_rivals)
Three clean blocks. Items 1 to 6 are the app and the website, items 8 to 14 are dealing with a person, items 15 to 18 are price and value. Those are defensible names: Digital experience, Support quality, Value for money.
Three items need a decision rather than a name.
q07, "it is easy to reach a person from inside the app", loads 0.55 and 0.53. That is not a measurement failure, it is the item doing exactly what it says: it is about the app and about support. Assigning it to one block and reporting the block as a clean measure of that construct would be a small lie. Leave it out of both indices and report it on its own.
q19, "compared with alternatives, this is good value", loads 0.62 on value and 0.43 on digital. Mostly value, but not cleanly, and the same treatment applies.
q20, "I enjoy the company's advertising", has a communality of 0.08. Ninety-two percent of its variance is unrelated to anything else in the battery. It is not measuring a fourth thing, it is measuring nothing this battery is about. Drop it from the indices and, if nobody uses it, from the questionnaire.
Cutting an item because it does not fit is a decision that should be made openly and written down, because it is also how a battery gets quietly tuned until it says what somebody wanted.
Step 8 · Reliability, and what alpha does not tell you¶
BLOCK = {"digital": Q[0:6], "support": Q[7:14], "value": Q[14:19]} # q07, q19, q20 excluded
def alpha(X):
X = np.asarray(X, float); k = X.shape[1]
return k/(k-1) * (1 - X.var(0, ddof=1).sum()/X.sum(1).var(ddof=1))
print("CRONBACH'S ALPHA")
for nm, cols in BLOCK.items():
print(f" {nm:9s} ({len(cols)} items) {alpha(e[cols]):.3f}")
print(f"\n all twenty items treated as one scale: {alpha(e[Q]):.3f}")
CRONBACH'S ALPHA digital (6 items) 0.839 support (7 items) 0.866 value (5 items) 0.828 all twenty items treated as one scale: 0.877
All three blocks are comfortably above the conventional 0.80, so each index is internally consistent.
Now look at the last line. Treating the whole battery as one scale gives an alpha of 0.877, higher than two of the three blocks, and we have just established that the battery measures three distinguishable things.
Alpha is not a test of unidimensionality. It rises with the number of items almost regardless of structure, so a high alpha on a long battery is close to uninformative. It answers "do these items hang together enough to average" and not "are these items measuring one thing".
Step 9 · What forgetting to reverse-score actually costs¶
Everything above used the reverse-scored data. Run the component analysis on the raw responses and compare.
p0, L0, _ = loadings(d[Q].astype(float))
print("eigenvalues, raw responses :", np.round(p0.explained_variance_[:5], 3))
print("eigenvalues, reverse-scored :", np.round(ev[:5], 3))
print("identical:", np.allclose(p0.explained_variance_, ev))
eigenvalues, raw responses : [6.174 2.559 2.178 0.952 0.659] eigenvalues, reverse-scored : [6.174 2.559 2.178 0.952 0.659] identical: True
Identical, and necessarily so. Reversing an item multiplies that column by minus one, which flips the sign of its correlations with everything and leaves the eigenvalues of the correlation matrix untouched.
That is exactly what makes the mistake dangerous. It does not announce itself. The scree plot looks fine, the parallel analysis gives the same answer, and the rotated solution still has three blocks. The only visible symptom is that four items carry negative loadings, which is easy to read past.
The damage is downstream, in the indices.
print("CRONBACH'S ALPHA, with and without reverse scoring")
for nm, cols in BLOCK.items():
print(f" {nm:9s} reverse-scored {alpha(e[cols]):.3f} raw {alpha(d[cols]):.3f}")
print("\nCORRELATION OF EACH INDEX WITH THE SEPARATE OVERALL SATISFACTION ITEM")
for nm, cols in BLOCK.items():
print(f" {nm:9s} reverse-scored {pearsonr(e[cols].mean(1), e.overall_satisfaction)[0]:+.3f}"
f" raw {pearsonr(d[cols].mean(1), e.overall_satisfaction)[0]:+.3f}")
CRONBACH'S ALPHA, with and without reverse scoring digital reverse-scored 0.839 raw 0.555 support reverse-scored 0.866 raw 0.246 value reverse-scored 0.828 raw 0.400 CORRELATION OF EACH INDEX WITH THE SEPARATE OVERALL SATISFACTION ITEM digital reverse-scored +0.545 raw +0.488 support reverse-scored +0.587 raw +0.470 value reverse-scored +0.532 raw +0.489
Alpha for the support block falls from 0.866 to 0.246. Within an index, a negatively worded item that has not been reversed points the opposite way to its neighbors and cancels them out.
So the practical guidance is not "remember to reverse-score", which everyone already knows. It is: compute alpha before you ship an index, because it is the diagnostic that catches this and the component analysis is not.
Step 10 · Do the weights earn their keep?¶
The dashboard needs a number it can recompute every month, in the reporting tool, without this notebook. Two candidates: the component scores, which need the loadings and the item means and standard deviations, or the plain unweighted average of the items in each block.
Zc = (e[Q].astype(float) - e[Q].astype(float).mean())/e[Q].astype(float).std(ddof=0)
scores = Zc.values @ LV @ np.linalg.inv(LV.T @ LV)
names = ["digital", "support", "value"]
print("COMPONENT SCORE AGAINST THE PLAIN AVERAGE OF THE SAME ITEMS")
for j, nm in enumerate(names):
print(f" {nm:9s} r = {pearsonr(scores[:,j], e[BLOCK[nm]].mean(1))[0]:+.4f}")
COMPONENT SCORE AGAINST THE PLAIN AVERAGE OF THE SAME ITEMS digital r = +0.9722 support r = +0.9769 value r = +0.9613
Above 0.96 for all three. The weights the analysis worked so hard to estimate are producing almost exactly the same ranking of respondents as adding the items up and dividing.
That is not a coincidence and it is not specific to this dataset. When items are positively correlated and their loadings are of similar size, almost any reasonable set of weights produces nearly the same composite. It is a result with a long history in decision research, usually filed under unit weights.
The question that actually settles it is whether either version predicts anything better.
kf = KFold(5, shuffle=True, random_state=2)
SUM = np.column_stack([e[c].mean(1) for c in BLOCK.values()])
cands = {"3 component scores": scores,
"3 unweighted averages": SUM,
"all 20 items separately": Zc.values}
print("five-fold cross-validation")
print(f" {'predictors':26s} {'overall satisfaction R2':>24s} {'would recommend AUC':>22s}")
for nm, X in cands.items():
r2 = cross_val_score(LinearRegression(), X, e.overall_satisfaction, cv=kf, scoring="r2").mean()
auc = cross_val_score(LogisticRegression(max_iter=2000), X, e.would_recommend, cv=kf,
scoring="roc_auc").mean()
print(f" {nm:26s} {r2:24.4f} {auc:22.4f}")
five-fold cross-validation predictors overall satisfaction R2 would recommend AUC 3 component scores 0.5696 0.7609 3 unweighted averages 0.5610 0.7605 all 20 items separately 0.5656 0.7552
Three numbers within 0.009 of each other on R-squared and within 0.006 on AUC.
Two conclusions follow, and they point in opposite directions about what the analysis was for.
Ship the averages. They are recomputable by anyone, they survive an item being added or dropped, they do not need to be re-estimated when next quarter's data arrives, and they cost nothing in predictive power. A component score that has to be regenerated by a data scientist every month is a worse product than a column of arithmetic in the reporting tool.
The analysis was still necessary. It is what established that there are three blocks and not one, which items belong to which, that q20 measures nothing relevant, and that q07 and q19 belong to two blocks at once. None of that was visible before, and averaging all twenty items into the single dashboard number marketing asked for would have hidden exactly the movement they wanted to see.
The dimension reduction was worth doing. The dimension reduction was not worth deploying.
Step 11 · The pictures¶
fig, axes = plt.subplots(1, 2, figsize=(12.6, 4.4))
ax = axes[0]
xs = np.arange(1, 11)
ax.plot(xs, ev[:10], color=VI, lw=2.8, marker="o", ms=6, label="this survey")
ax.plot(xs, p95[:10], color=RD, lw=2.0, ls="--", marker="s", ms=4, label="95th percentile of noise")
ax.axhline(1, color=MUT, lw=1.6, ls=":", label="Kaiser cutoff of 1")
ax.set_xticks(xs); ax.set_xlabel("component"); ax.set_ylabel("eigenvalue")
ax.set_title("Three beat noise, and the fourth does not")
ax.legend(fontsize=9)
ax = axes[1]
cols = [VI, LT, MUT]
for j, (nm, cs) in enumerate(BLOCK.items()):
s_ = e[cs].mean(1).values
r = np.corrcoef(scores[:, j], s_)[0, 1]
ax.scatter(s_, scores[:, j], s=7, alpha=0.30, color=cols[j], label=f"{nm} r = {r:.3f}")
ax.set_xlabel("plain average of the block")
ax.set_ylabel("component score")
ax.set_title("The weights and the arithmetic agree")
leg = ax.legend(fontsize=9, loc="upper left", markerscale=2.4)
for h in leg.legend_handles: h.set_alpha(1)
plt.tight_layout(); plt.show()
Left: the retention decision, with the noise reference that makes the eigenvalues interpretable. Right: the component scores against the plain averages we are recommending instead of them.
Step 12 · Two things this does not settle¶
PCA is not factor analysis, and the difference matters if the goal changes. Principal components are weighted sums of the observed items, chosen to capture variance. A factor model runs the other way: it treats each item as a noisy indicator of an unobserved construct plus item-specific error, and estimates the construct. For building indices, which is the job here, PCA is the right tool and the two give near-identical answers when communalities are high. For a claim like "customer trust causes retention", where trust is a latent thing the items only indicate, the factor model is the object that means what you want it to mean.
Pearson correlations on five-point items understate the association. Everything above treats a Likert response as a number, which is standard and slightly wrong: the responses are ordered categories, and discretising an underlying continuum attenuates correlations, more so as the items get more skewed. Polychoric correlations estimate the correlation of the latent continuous variables behind the categories and would give slightly larger loadings here. With five points and mild skew the difference rarely changes which items group together, which is why the shortcut is usually taken. With three-point items or heavy floor effects, it can.
What to take away¶
- Kaiser's rule retains about ten components from twenty variables with no structure at all. The threshold of 1 is the average eigenvalue, not a test.
- Parallel analysis compares against what noise produces at this sample size, and it is the rule to reach for by default. It rejected the fourth component by 0.18 where Kaiser was within 0.05 of keeping it.
- Rotation does not change the fit and does not find anything. Same 54.6 percent of variance, different axes, and only the rotated ones are readable.
- Not every item belongs somewhere. Two cross-load and one has a communality of 0.08; the honest move is to say so rather than to force them into a block.
- A high alpha is not evidence of one construct. All twenty items together score 0.877, and they are three things.
- The plain average of each block correlates above 0.96 with the component score and predicts both outcomes as well. Do the analysis; ship the arithmetic.