9 Multivariate Statistical Methods

Many of the things researchers most want to study cannot be measured directly. Stress, burnout, satisfaction, intelligence, motivation, and quality of life are constructs (Chapter 5): ideas that exist only through their effects on what people say and do. A questionnaire approaches such a construct indirectly, through several items that are each expected to reflect it, and the researcher then has to show that the items really do measure what they are meant to. At the same time, a dataset with many variables often contains patterns that no single variable reveals, such as groups of people with a similar profile across all of them.

Multivariate methods analyse many variables at once, and this chapter introduces three of them. Principal component analysis summarises many variables with a few. Factor analysis checks whether questionnaire items measure the underlying traits they are meant to, and Cronbach’s alpha checks whether they do so consistently. Cluster analysis sorts people into groups with similar profiles. In the study, they answer two research questions: whether the 22 questionnaire items measure stress, burnout, supervisor support, and satisfaction as intended (RQ6), a question every examiner of a questionnaire study will ask, and whether there are distinct profiles of students (RQ7).

TipBy the end of this chapter you will be able to
  • Explain why a construct is measured with several items, and what a latent variable is.
  • Explore the correlations among many variables at once.
  • Run a principal component analysis in Python, and decide how many components to keep.
  • Explain the difference between principal component analysis and factor analysis.
  • Run an exploratory factor analysis, and interpret factor loadings, cross-loadings, and reversed items.
  • Check the reliability of a questionnaire scale with Cronbach’s alpha.
  • Group cases with k-means and hierarchical clustering, choose the number of clusters, and describe the clusters as summaries rather than natural groups.

9.1 Constructs and latent variables

A single questionnaire item is a poor measure of a construct. The answer to “I feel unable to control important things in my studies” depends on the student’s stress, but also on how they read the question that day, on their mood, and on how they use the answer scale. Each item therefore carries two things: a signal from the construct, and noise of its own. A latent variable is the construct thought to lie behind the items: it is not observed, but it is assumed to cause part of every answer.

Averaging several items keeps the signal, which all items share, and lets the noise, which differs from item to item, partly cancel out. A simulation shows how much this helps. It creates a “true” stress level for 600 imaginary students, then six items, each equal to the true level plus its own random noise, and compares how closely one item and the average of several items follow the true level:

np.random.seed(42)
true_stress = np.random.normal(0, 1, 600)
simulated_items = np.column_stack([true_stress + np.random.normal(0, 1, 600) for _ in range(6)])

corrs = {}
for k in range(1, 7):
    mean_k = simulated_items[:, :k].mean(axis=1)
    r = float(np.corrcoef(mean_k, true_stress)[0, 1])
    corrs[f"{k} items"] = round(r, 2)

for label, val in corrs.items():
    print(f"{label}: r = {val}")
1 items: r = 0.68
2 items: r = 0.8
3 items: r = 0.85
4 items: r = 0.88
5 items: r = 0.9
6 items: r = 0.92

In this simulation, one item correlates only about 0.72 with the true stress level, while the average of six items correlates about 0.93. This is why questionnaires use several items for each construct, and why the items are combined into a scale score (Chapter 3).

The approach rests on an assumption that must be checked: that the items of a scale really reflect one common construct, and not several. If some stress items actually measured tiredness, averaging them would mix two things. Chapter 5 called this construct validity. The methods in this chapter provide evidence for it: factor analysis shows whether the items group as intended, and Cronbach’s alpha shows whether the items of each group agree with each other, which Chapter 5 called internal consistency.

9.2 Many variables at once

The 22 questionnaire items are in questionnaire:

items = questionnaire.drop(columns=["student_id"])
print(f"Number of questionnaire items: {items.shape[1]}")
Number of questionnaire items: 22

With 22 items, there are 231 correlations between pairs of items, far too many to read one by one. A picture helps. In Python, we can sort the correlation matrix using hierarchical clustering so that items that correlate strongly sit together:

corr = items.corr()

# Sort items using hierarchical clustering order
dists = squareform(1 - corr.to_numpy())
order = leaves_list(linkage(dists, method="average"))
sorted_corr = corr.iloc[order, order]

fig, ax = plt.subplots(figsize=(7, 6.5))
sns.heatmap(
    sorted_corr,
    cmap="vlag",
    vmin=-1,
    vmax=1,
    center=0,
    square=True,
    cbar_kws={"shrink": 0.75, "label": "Correlation (r)"},
    ax=ax
)
plt.xticks(rotation=90, fontsize=8)
plt.yticks(rotation=0, fontsize=8)
plt.tight_layout()
plt.show()
A 22 by 22 grid of coloured squares. Four blocks of strongly correlated items stand out along the diagonal, one for each scale; stress and burnout items also correlate with each other.
Figure 9.1: Correlations among the 22 questionnaire items, sorted so that similar items are next to each other.

Blocks of related items stand out along the diagonal, one for each scale, just as the idea of latent variables predicts: items that share a construct correlate with each other. Blue squares show positive correlations and red negative. The stress and burnout blocks also correlate with each other, and one item, stress_4, is red against the other stress items, because it is worded the other way round. The methods in this chapter turn this picture into numbers.

9.3 Principal component analysis

Principal component analysis (PCA) replaces many correlated variables with a few new ones, called principal components, that together keep most of the information. The first component is the combination of the variables that captures as much of their variation as possible. The second captures as much as possible of what is left, and so on.

Two items that correlate strongly, such as “I feel emotionally drained by my studies” and “I feel exhausted when I think about my thesis”, illustrate the idea. Most students who score high on one score high on the other, so one combined score, their average for example, keeps most of what the two items tell. PCA does this for all the variables at once, and finds the best combinations automatically.

9.3.1 Running PCA

In Python, scikit-learn provides PCA, and StandardScaler standardises every variable to unit variance first, which is standard practice in survey analysis. Like classic PCA algorithms, standard PCA needs complete data, so .dropna() keeps students who answered every item:

items_complete = items.dropna()
print(f"Complete cases: {len(items_complete)} out of {len(items)}")

scaler = StandardScaler()
items_scaled = scaler.fit_transform(items_complete)

pca = PCA()
pca.fit(items_scaled)

# Variance explained by each component
exp_var = pca.explained_variance_ratio_
cum_var = np.cumsum(exp_var)
eigenvalues = pca.explained_variance_

pca_summary = pd.DataFrame({
    "Component": [f"PC{i+1}" for i in range(len(exp_var))],
    "Eigenvalue": eigenvalues.round(3),
    "Proportion": exp_var.round(3),
    "Cumulative": cum_var.round(3),
})
pca_summary.head(6)
Complete cases: 396 out of 600
Component Eigenvalue Proportion Cumulative
0 PC1 6.672 0.302 0.302
1 PC2 2.876 0.130 0.433
2 PC3 1.567 0.071 0.504
3 PC4 1.234 0.056 0.560
4 PC5 0.809 0.037 0.597
5 PC6 0.727 0.033 0.630

Only 396 of the 600 students answered all 22 items: a small share of missing answers per item adds up to many incomplete students. Factor analysis, below, can use the incomplete students too.

The table shows how much of the total variation each component captures. The first component alone captures 30%, and the first four together capture 56%.

9.3.2 The number of components

No single rule decides how many components to keep, so researchers use several. The Kaiser rule keeps components whose eigenvalue, the amount of variation they capture, is greater than 1, that is, more than one original variable’s worth; here, 4 components pass. The scree plot shows the eigenvalues in order, and the researcher looks for the “elbow” where the line flattens out. Parallel analysis compares the eigenvalues with those from random data of the same size, and keeps the components that beat random data; it is the most reliable of the three.

We can draw the scree plot directly with Matplotlib and Seaborn:

fig, ax = plt.subplots(figsize=(6.5, 3.8))
components_10 = range(1, 11)
ax.bar(components_10, exp_var[:10] * 100, color=navy, alpha=0.85, label="% of Variance")
ax.plot(components_10, exp_var[:10] * 100, color=copper, marker="o", linewidth=2.2, label="Trend line")

for i, pct in enumerate(exp_var[:10] * 100, 1):
    ax.text(i, pct + 0.6, f"{pct:.1f}%", ha="center", fontsize=8)

ax.set_xlabel("Principal Component")
ax.set_ylabel("Variance Explained (%)")
ax.set_xticks(components_10)
ax.set_xticklabels([f"PC{i}" for i in components_10])
ax.set_ylim(0, max(exp_var[:10] * 100) + 4)
plt.tight_layout()
plt.show()
Bar chart with a line, falling steeply from the first component to the fourth, then flattening from the fifth onwards.
Figure 9.2: Scree plot: the share of variation captured by each principal component.

The line drops steeply for four components and flattens after that. Parallel analysis can be run by simulating random Gaussian datasets of the same dimensions:

def parallel_analysis(data, n_iter=100, seed=42):
    np.random.seed(seed)
    n, p = data.shape
    rand_eigs = np.zeros((n_iter, p))
    for i in range(n_iter):
        synth = np.random.normal(0, 1, size=(n, p))
        pca_rand = PCA().fit(StandardScaler().fit_transform(synth))
        rand_eigs[i] = pca_rand.explained_variance_
    return rand_eigs.mean(axis=0)

rand_eigs = parallel_analysis(items_scaled)
n_factors = int(np.sum(eigenvalues > rand_eigs))
print(f"Parallel analysis suggests {n_factors} components.")
Parallel analysis suggests 3 components.

All three methods point to four dimensions, matching the four scales the questionnaire was designed to measure.

9.3.3 PCA and factor analysis compared

PCA and factor analysis are often confused, and many theses use one when they mean the other. The difference lies in what they assume. PCA simply summarises the variables: the components are combinations of the items, with no claim about why the items are related, and it is used to reduce many variables to a few, for example before a further analysis. Factor analysis assumes the model of the first section: that each item is caused by one or more latent variables, called factors, plus error of its own. A student answers the stress items the way they do because they are stressed. Factor analysis is therefore the method for checking whether a questionnaire measures the constructs it is meant to, and that is the study’s question.

9.4 Factor analysis

Exploratory factor analysis (EFA) estimates how strongly each item is related to each factor. These relationships are called loadings. Loadings run roughly from −1 to 1: close to 0 means that the item is unrelated to the factor, and a value above about 0.4 in size is usually considered a clear relationship.

In Python, statsmodels.multivariate.factor.Factor provides factor analysis with multiple estimation techniques and rotations. Promax or oblimin rotation allows the factors to correlate with each other, which is realistic, since stress and burnout surely go together:

# Fit exploratory factor analysis with 4 factors
fa = Factor(items_complete, n_factor=4, method="pa").fit()
fa.rotate("promax")

loadings = pd.DataFrame(
    fa.loadings,
    index=items_complete.columns,
    columns=["Factor 1", "Factor 2", "Factor 3", "Factor 4"]
)

# Display loadings for the first 10 items
loadings.head(10).round(2)
Factor 1 Factor 2 Factor 3 Factor 4
stress_1 -0.65 0.19 -0.12 -0.23
stress_2 -0.51 0.20 -0.09 -0.33
stress_3 -0.55 0.09 -0.07 -0.33
stress_4 0.52 -0.11 0.08 0.27
stress_5 -0.69 0.13 -0.17 -0.26
stress_6 -0.51 0.08 -0.12 -0.21
burnout_1 -0.31 0.11 -0.14 -0.66
burnout_2 -0.21 0.00 -0.02 -0.64
burnout_3 -0.46 0.14 -0.06 -0.56
burnout_4 -0.17 0.07 -0.15 -0.63

The table shows four clear factors. Each collects the items of one scale: the support items on one, the stress items on another, the burnout items on a third, and the satisfaction items on the fourth. The reversed item, stress_4 (“I feel confident handling problems in my studies”), loads with the opposite sign (+0.52 against −0.51 to −0.69 for the other stress items), which is exactly right for a reversed item: agreeing with it means less stress, and it is why the item must be reversed before the stress score is calculated (Chapter 3). One item, burnout_3 (“Deadlines make me feel overwhelmed”), loads substantially on both burnout (−0.56) and stress (−0.46). Reading the item, that makes sense: feeling overwhelmed by deadlines is both. Such a cross-loading item is worth discussing in a thesis, and sometimes worth rewording or removing in future studies.

9.4.1 Reliability and Cronbach’s alpha

Factor analysis shows that the items of each scale belong together. Reliability asks a related question: whether the items of a scale give consistent results (Chapter 5). The most common measure of this internal consistency is Cronbach’s alpha, which ranges from 0 to 1:

\[ \alpha = \frac{k}{k - 1} \left(1 - \frac{\sum s_i^2}{s_{\text{total}}^2}\right) \]

It rises when the items correlate strongly with each other, and also when there are more items, for the reason shown by the simulation at the start of the chapter. As a rough guide, 0.7 or above is acceptable and 0.8 or above good for research use:

def cronbach_alpha(df: pd.DataFrame) -> float:
    """Calculate Cronbach's alpha for internal consistency."""
    df_clean = df.dropna()
    k = df_clean.shape[1]
    item_variances = df_clean.var(ddof=1).sum()
    total_variance = df_clean.sum(axis=1).var(ddof=1)
    return float((k / (k - 1)) * (1 - (item_variances / total_variance)))

q = questionnaire.copy()
q["stress_4"] = 6 - q["stress_4"]  # reverse item 4

alphas = {
    "Stress": cronbach_alpha(q.filter(like="stress_")),
    "Burnout": cronbach_alpha(q.filter(like="burnout_")),
    "Support": cronbach_alpha(q.filter(like="support_")),
    "Satisfaction": cronbach_alpha(q.filter(like="satisfaction_")),
}

for scale_name, alpha_val in alphas.items():
    print(f"{scale_name:12s}: alpha = {alpha_val:.2f}")
Stress      : alpha = 0.83
Burnout     : alpha = 0.83
Support     : alpha = 0.85
Satisfaction: alpha = 0.76

All four scales are reliable, with alphas between 0.76 and 0.85.

A high alpha shows consistency, not validity. A scale can be highly reliable and still measure the wrong thing, as the target picture in Chapter 5 showed; the factor structure and the relationships with other constructs are the evidence for validity.

TipReporting the questionnaire in a thesis

A typical report: “An exploratory factor analysis (principal axis factoring, promax rotation) supported four factors, matching parallel analysis and scree plot criteria. All items loaded predominantly on their intended factor (loadings > 0.40), with one item (burnout_3) cross-loading on both burnout and stress. Internal consistency was acceptable to good across all scales (Cronbach’s \(\alpha = 0.76\) to \(0.85\)).”

9.5 Cluster analysis

Factor analysis groups variables. Cluster analysis groups people: it looks for groups of students whose profiles are similar to each other and different from other groups. It is a descriptive method. It does not test a hypothesis, and the groups it finds are summaries of the data, useful for describing and thinking about students, not discoveries of natural kinds of student. The study’s question about student profiles (RQ7) is accordingly exploratory (Chapter 5).

Each student is described by their average sleep, study hours, caffeine, and exercise across the semesters, and their stress, support, and satisfaction scores:

scores = pd.DataFrame({
    "student_id": q["student_id"],
    "stress": q.filter(like="stress_").mean(axis=1),
    "support": q.filter(like="support_").mean(axis=1),
    "satisfaction": q.filter(like="satisfaction_").mean(axis=1),
})

profiles = (
    semesters.groupby("student_id")[["sleep_hours", "study_hours", "caffeine_mg", "exercise_days"]]
    .mean()
    .reset_index()
    .merge(scores, on="student_id")
    .dropna()
)

scaler = StandardScaler()
profile_features = profiles.drop(columns=["student_id"])
profile_scaled = scaler.fit_transform(profile_features)

Clustering works with distances between students: how different two students’ profiles are. Caffeine is measured in hundreds of milligrams and sleep in single hours, so without adjustment caffeine would dominate every distance. StandardScaler turns every variable into z-scores (Chapter 6), so each counts equally.

9.5.1 k-means clustering

k-means is the most widely used clustering method. The researcher chooses the number of clusters, \(k\), and the algorithm places \(k\) starting points at random, assigns each student to the nearest point, moves each point to the centre (the mean) of its students, and repeats the last two steps until nothing changes. In scikit-learn, n_init=25 runs the algorithm from 25 different starting seeds and keeps the best partition.

9.5.2 The number of clusters

As with components, no single rule decides the number of clusters, and two guides are common. The elbow method uses the fact that the total distance of students from their cluster centres (inertia) falls as \(k\) grows, and looks for the point where adding clusters stops helping much. The silhouette measures, for each student, how much closer they are to their own cluster than to the next nearest one, from −1 to 1; a higher average silhouette means clearer clusters.

inertias = []
silhouettes = []
k_range = range(2, 9)

for k in k_range:
    km = KMeans(n_clusters=k, n_init=25, random_state=123).fit(profile_scaled)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(profile_scaled, km.labels_))

fig, axes = plt.subplots(1, 2, figsize=(8, 3.4))

# Elbow plot
axes[0].plot(k_range, inertias, marker="o", color=navy, linewidth=2)
axes[0].set_xlabel("Number of clusters (k)")
axes[0].set_ylabel("Inertia (Total within-cluster sum of squares)")
axes[0].set_title("Elbow Method")

# Silhouette plot
axes[1].plot(k_range, silhouettes, marker="s", color=copper, linewidth=2)
axes[1].set_xlabel("Number of clusters (k)")
axes[1].set_ylabel("Average Silhouette Score")
axes[1].set_title("Silhouette Method")

plt.tight_layout()
plt.show()
Two plots. Left: a falling line of within-cluster inertia. Right: average silhouette width peaking at 2 clusters and remaining moderate for 3 and 4.
Figure 9.3: Choosing the number of clusters: inertia elbow method (left) and silhouette score (right).

The honest conclusion is that the data does not point clearly to one number. The elbow is gentle, and the silhouette is highest for 2 clusters but only a little lower for 3 or 4. This is common with real data: groups of people overlap, and there are no sharp boundaries. The choice then rests on what is useful and interpretable, and on prior theory. The thesis, based on earlier research, expects around four student profiles, so \(k = 4\) is tried:

km4 = KMeans(n_clusters=4, n_init=25, random_state=123)
profiles["cluster"] = km4.fit_predict(profile_scaled) + 1
profiles["cluster"].value_counts().sort_index()
cluster
1     90
2    201
3    113
4    195
Name: count, dtype: int64

9.5.3 Describing the clusters

A cluster is only useful if its meaning can be stated. Adding the cluster to the original, unscaled data and averaging each variable by cluster shows each group’s profile:

cluster_summary = (
    profiles.groupby("cluster")
    .agg(
        students=("student_id", "count"),
        sleep_hours=("sleep_hours", "mean"),
        study_hours=("study_hours", "mean"),
        caffeine_mg=("caffeine_mg", "mean"),
        exercise_days=("exercise_days", "mean"),
        stress=("stress", "mean"),
        support=("support", "mean"),
        satisfaction=("satisfaction", "mean")
    )
    .round(1)
    .reset_index()
)
cluster_summary
cluster students sleep_hours study_hours caffeine_mg exercise_days stress support satisfaction
0 1 90 5.2 40.9 414.0 1.3 3.7 2.9 2.9
1 2 201 7.2 23.7 123.2 3.5 2.6 3.6 3.8
2 3 113 5.8 41.2 187.9 1.5 3.6 3.5 3.3
3 4 195 6.7 17.9 148.0 2.2 3.4 2.7 2.5

The averages describe the groups. One is clearly balanced: the most sleep, moderate study, the most exercise, low stress, and good support. Two are overloaded: long study weeks, little sleep and exercise, high stress, and one of them with very high caffeine intake. The fourth combines low study hours with low support and low satisfaction: students who seem disengaged or isolated. Cluster numbers are arbitrary labels, and can change if the code is run with a different seed.

Figure 9.4 shows the clusters on the first two principal components of the profile data, a common way to see many variables in two dimensions:

pca_2d = PCA(n_components=2).fit_transform(profile_scaled)
profiles["PC1"] = pca_2d[:, 0]
profiles["PC2"] = pca_2d[:, 1]

fig, ax = plt.subplots(figsize=(7, 4.5))
sns.scatterplot(
    data=profiles,
    x="PC1",
    y="PC2",
    hue="cluster",
    palette="Set2",
    alpha=0.75,
    s=45,
    ax=ax
)
ax.set_xlabel("First Principal Component (PC1)")
ax.set_ylabel("Second Principal Component (PC2)")
ax.legend(title="Cluster", bbox_to_anchor=(1.02, 1), loc="upper left")
plt.tight_layout()
plt.show()
Scatter plot of students on two principal components, coloured by cluster, with four overlapping groups.
Figure 9.4: The four k-means clusters, shown on the first two principal components of the profile data.

The clusters overlap. That is not a failure: real students do not fall into neat boxes. The profiles are useful summaries, not natural categories.

9.5.4 Hierarchical clustering

Hierarchical clustering takes a different approach. It starts with every student as their own cluster and repeatedly merges the two most similar clusters, until everyone is in one. The result is a tree, called a dendrogram, which can be cut at any height to give any number of clusters, without choosing \(k\) in advance. Ward’s method merges clusters so that they stay as compact as possible:

Z = linkage(profile_scaled, method="ward")

fig, ax = plt.subplots(figsize=(7.5, 4))
dendrogram(
    Z,
    truncate_mode="lastp",
    p=25,
    show_leaf_counts=True,
    color_threshold=Z[-3, 2],
    ax=ax
)
ax.set_xlabel("Cluster branches (truncated)")
ax.set_ylabel("Ward Distance")
plt.tight_layout()
plt.show()
A tree diagram whose branches merge upwards; coloured branches mark clusters where the tree is cut.
Figure 9.5: Dendrogram from hierarchical clustering (Ward’s method), cut into four clusters.

We can cut the hierarchical tree into four clusters with fcluster and compare it with k-means:

hier_clusters = fcluster(Z, t=4, criterion="maxclust")
table_agreement = pd.crosstab(hier_clusters, profiles["cluster"], rownames=["Hierarchical"], colnames=["k-means"])
table_agreement
k-means 1 2 3 4
Hierarchical
1 8 10 74 1
2 74 0 33 1
3 8 4 3 150
4 0 187 3 43

They agree on the broad picture, most clearly on the balanced students, but many students switch groups between the two methods. When two reasonable methods disagree about a student, that student is near a boundary. The clearer groups, found by both methods, are the ones to trust most.

9.5.5 Clusters in random data

A clustering method will split any data into groups, even data that contains no groups at all. The final code makes the point with 300 points drawn entirely at random, with no structure of any kind, and asks k-means for four clusters:

np.random.seed(7)
random_points = pd.DataFrame(np.random.normal(0, 1, (300, 2)), columns=["x", "y"])
km_rand = KMeans(n_clusters=4, n_init=25, random_state=7).fit(random_points)
random_points["cluster"] = km_rand.labels_ + 1

fig, ax = plt.subplots(figsize=(5, 4))
sns.scatterplot(
    data=random_points,
    x="x",
    y="y",
    hue="cluster",
    palette="viridis",
    s=40,
    legend=False,
    ax=ax
)
ax.set_aspect("equal")
ax.set_title("Clusters Found in Pure Noise")
plt.tight_layout()
plt.show()
Scatter plot of 300 points in a single round cloud, coloured in four wedge-shaped regions by the cluster k-means assigned.
Figure 9.6: k-means finds four ‘clusters’ in 300 random points that contain no groups at all.

The method obliges, and divides one round cloud into four neat regions. Nothing in its output says that the groups are artificial. The existence of clusters therefore proves nothing on its own. Before clusters are reported as meaningful, it should be checked that they make sense, that different methods broadly agree, and that the groups differ in ways that matter for the research. Chapter 14 introduces methods that allow for overlapping groups and for students who fit no group at all.

NoteIn your field: criminology and social policy

The classic USArrests dataset records rates of assault, murder, and rape per 100,000 residents, and the percentage of people living in urban areas, for the 50 US states in 1973. PCA summarises the four variables:

X_arrests = StandardScaler().fit_transform(usarrests.drop(columns=["State"]))
pca_arr = PCA().fit(X_arrests)
var_arr = pca_arr.explained_variance_ratio_

for i, v in enumerate(var_arr[:2], 1):
    print(f"PC{i}: {v * 100:.1f}% of variance")
PC1: 62.0% of variance
PC2: 24.7% of variance

The first component captures about 62% of the variation: it is essentially an overall violent crime rate. The second, about 25%, mostly reflects how urban a state is. Two numbers now describe what took four.

9.6 Common misconceptions

Multivariate methods produce impressive output, which makes their limits easy to forget.

  • “PCA and factor analysis are the same.” PCA summarises variables; factor analysis models latent variables that cause them. Only factor analysis addresses whether a questionnaire measures its constructs.
  • “A high Cronbach’s alpha proves a scale is valid.” Alpha measures consistency. A scale can be consistent and still measure the wrong thing.
  • “Clusters are natural groups.” Clustering finds groups in any data, including random data. Clusters are summaries, and their usefulness must be argued.
  • “The number of factors or clusters is decided by the software.” The criteria often disagree, and the final choice is a judgement that should be justified in the thesis.

9.7 Chapter review

9.7.1 Summary

  • A construct is measured with several items because each item carries noise of its own; averaging keeps the shared signal. A latent variable is the unobserved construct assumed to cause the items.
  • A sorted correlation plot shows the structure among many variables.
  • PCA replaces many correlated variables with a few components that keep most of the variation. Decide how many to keep with eigenvalues above 1, the scree plot, and, most reliably, parallel analysis.
  • PCA summarises; factor analysis models latent variables that cause the items. Use factor analysis to check a questionnaire.
  • In an EFA, loadings show how strongly items relate to factors. Check that items load on their intended factor, that reversed items load with opposite signs, and discuss cross-loadings. Oblique rotation lets factors correlate.
  • Cronbach’s alpha measures a scale’s internal consistency: 0.7 or above is acceptable, 0.8 or above good. Reverse reversed items first. Reliability is necessary for validity but not enough.
  • Cluster analysis groups people with similar profiles. Scale the variables first. k-means needs the number of clusters in advance; hierarchical clustering produces a tree that can be cut anywhere. Choose the number with the elbow and silhouette methods, theory, and interpretability.
  • Clustering finds groups even in random data. Clusters are useful summaries, not proof of natural groups; describe each cluster by its averages, and check that methods agree.

9.7.2 Key terms

Multivariate analysis, construct, latent variable, correlation matrix, principal component analysis, principal component, eigenvalue, Kaiser rule, scree plot, parallel analysis, factor analysis, factor, loading, rotation, oblique rotation, cross-loading, reversed item, reliability, internal consistency, Cronbach’s alpha, confirmatory factor analysis, cluster analysis, distance, scaling, k-means, elbow method, silhouette, hierarchical clustering, dendrogram, Ward’s method.

9.8 Exercises

The playground has these and more, with hints and solutions.

  1. Run a PCA on the six support items only. Report how much variation the first component captures, and what that suggests about the scale.
  2. Run an EFA with three factors instead of four. Identify the scales that end up sharing a factor, and suggest why.
  3. Calculate Cronbach’s alpha for the stress scale without reversing stress_4 first, and explain what happens.
  4. Run k-means with 3 clusters on the profile data, describe the clusters, and identify which profiles from the 4-cluster solution were merged.
  5. Explain, in two or three sentences for a thesis methods section, why you chose the number of clusters you did.
  6. Change the simulation at the start of the chapter so that each item has more noise (sd = 2). Find how many items are now needed for the average to correlate at least 0.8 with the true stress level.

9.9 Further reading

  • An Introduction to Statistical Learning: with Applications in Python (James et al. 2023), chapter “Unsupervised Learning”, covers PCA and clustering in depth with Python implementations.
  • Applied Multivariate Statistical Analysis (Johnson and Wichern 2007) is the classic comprehensive reference for factor analysis, PCA, and clustering theory.

References

James, Gareth, Daniela Witten, Trevor Hastie, Robert Tibshirani, and Jonathan Taylor. 2023. An Introduction to Statistical Learning: With Applications in Python. Springer. https://doi.org/10.1007/978-3-031-38747-0.
Johnson, Richard A., and Dean W. Wichern. 2007. Applied Multivariate Statistical Analysis. 6th ed. Pearson Prentice Hall.