Playground: Chapter 9
Multivariate Statistical Methods
This page practises the ideas of Chapter 9 in four parts. Practise the chapter works through the book’s exercises, with starter code, hints, and solutions. Go further goes beyond the book, with simulations and real datasets. Check your understanding asks short questions whose answers open with a click. Do it yourself gives open tasks with an empty code space and no answers: you decide how to check and summarise the data, as you will in your own thesis.
To keep the page fast, it uses standard Python data science tools: scikit-learn (PCA, StandardScaler, KMeans), statsmodels (Factor), and pandas/numpy for Cronbach’s alpha and summary matrices. In the exercises, replace each ______ with your own code and press Run Code.
Practise the chapter
These are the exercises at the end of Chapter 9, with the same numbers.
Exercise 1: PCA of the support items
Run a PCA on the six support items only. Report how much variation the first component captures, and what that suggests about the scale.
In scikit-learn, fitting a transformer or estimator to data uses .fit().
support_items = items.filter(like="support_").dropna()
scaler_sup = StandardScaler()
support_scaled = scaler_sup.fit_transform(support_items)
pca_support = PCA()
pca_support.fit(support_scaled)
print("Explained variance ratio:", np.round(pca_support.explained_variance_ratio_, 3))The first component captures 57% of the variation in the six items, and each of the other five between 7% and 10%. One strong component and five weak ones mean that the items mostly measure one thing, which is what a scale should do: averaging them into one support score loses little. (.dropna() kept the 537 students who answered all six items.)
Exercise 2: Three factors instead of four
Run a factor analysis with three factors instead of four. Identify the scales that end up sharing a factor, and suggest why.
n_factor specifies the number of factors to extract (here, 3).
items_complete = items.dropna()
fa3 = Factor(items_complete, n_factor=3, method="pa").fit()
fa3.rotate("promax")
loadings3 = pd.DataFrame(fa3.loadings, index=items.columns, columns=["F1", "F2", "F3"])
print(loadings3.mask(loadings3.abs() < 0.3, "").to_string())All six stress items and all six burnout items load on the first factor (stress_4 negatively, since it is not reversed here); support and satisfaction keep a factor each. Stress and burnout are closely related, since students who feel overwhelmed also feel worn out, and their factors correlated strongly in the chapter, so when only three factors are allowed, they are the pair that merges. Parallel analysis in the chapter supported four factors, and the four-factor solution keeps them apart. Note that dropping missing rows leaves the 396 students who answered every item.
Exercise 3: Cronbach’s alpha without reversing
Calculate Cronbach’s alpha for the stress scale without reversing stress_4 first, and explain what happens. The function below calculates alpha from its formula: the number of items, \(k\), the variance of each item, and the variance of the total score:
\[ \alpha = \frac{k}{k - 1}\left(1 - \frac{\sum \text{item variances}}{\text{variance of the total}}\right) \]
On a 1-to-5 Likert scale, subtracting an answer from 6 reverses it (1 becomes 5, 5 becomes 1).
stress_reversed = stress_items.copy()
stress_reversed["stress_4"] = 6 - stress_reversed["stress_4"]
print("Reversed alpha:", round(cronbach_alpha(stress_reversed), 2))Without reversing, alpha is 0.57, below the conventional threshold of 0.70; with stress_4 reversed, it is 0.83, good. The formula shows why: alpha is high when the total varies much more than the items separately, which happens when the items rise and fall together. An unreversed item moves against the others, cancels part of their shared variation, and makes the total less variable.
Exercise 4: Three clusters
Run k-means with 3 clusters on the standardized profile data, describe the clusters, and identify which profiles from the 4-cluster solution were merged.
n_clusters sets the target number of clusters (here, 3). In the cross-tabulation table, a column of the four-cluster solution whose students all fall into one row has been absorbed into that three-cluster group.
kmeans3 = KMeans(n_clusters=3, n_init=25, random_state=123)
profiles_k3 = profiles.copy()
profiles_k3["cluster_3"] = kmeans3.fit_predict(profile_scaled)
print(profiles_k3.groupby("cluster_3").mean().round(1))
kmeans4 = KMeans(n_clusters=4, n_init=25, random_state=123)
profiles_k3["cluster_4"] = kmeans4.fit_predict(profile_scaled)
print(pd.crosstab(profiles_k3["cluster_3"], profiles_k3["cluster_4"]))The three clusters are the balanced students (about 208: 7.2 hours of sleep, low stress, high support and satisfaction), the disengaged (about 206: only 19 hours of study a week, low support and satisfaction), and the overloaded (about 185: 42 hours of study, 5.5 hours of sleep, and a high average caffeine intake). The cross-tabulation shows that the overloaded cluster absorbed two clusters of the four-cluster solution: overloaded students with very high caffeine (over 400 mg a day) and those with moderate caffeine. With three clusters, caffeine no longer divides the overloaded students.
Exercise 5: Justifying the number of clusters
Explain, in two or three sentences for a thesis methods section, why you chose the number of clusters you did. Write your answer first, then open the model answer.
“The elbow and silhouette criteria did not identify an unambiguous number of clusters: the average silhouette was highest for two clusters but only slightly lower for three and four. We therefore chose four clusters, the number of student profiles expected from previous higher-education research, because the four clusters were interpretable and each contained a substantial share of the sample; a three-cluster solution merged two overloaded profiles that differed mainly in caffeine intake. The clusters are reported as descriptive summaries, not as natural types.” The paragraph names the criteria, admits that they were not decisive, gives the substantive theoretical reason for the choice, and clarifies how the clusters should be interpreted.
Exercise 6: Noisier items
Change the simulation at the start of the chapter so that each item has more noise (scale=2 instead of 1), and use up to 20 items. Find how many items are now needed for the average to correlate at least 0.8 with the true stress level.
The noise of each item is set by the scale argument of np.random.normal(); the chapter used 1.
simulated_items = np.column_stack([true_stress + np.random.normal(scale=2, size=600) for _ in range(20)])With twice the noise, one item correlates only ~0.45 with the true stress level, and the average needs 8 items to reach 0.80; six items, enough in the chapter, now reach only ~0.76. Theory agrees: with noise twice as large as the signal, the correlation of the average of \(k\) items is \(1/\sqrt{1 + 4/k}\), which passes 0.8 at \(k = 8\). Noisy items can be rescued by using more of them, but each extra item helps less than the one before.
Go further
These exercises go beyond the book.
Exercise 7: When alpha misleads
Alpha measures how much items share, not whether what they share is the construct. Simulate six stress items as in the chapter, then add the same extra noise to every item, a “mood” that colours every answer on the day of the survey. Compare alpha, and the correlation of the average with the true stress level, with and without the mood.
Adding mood[:, np.newaxis] broadcasts the 600-element vector across all 6 columns of the simulated matrix.
items_mood = items_clean + mood[:, np.newaxis]With the mood, alpha rises, from ~0.86 to ~0.92, while the average’s correlation with the true stress level falls, from ~0.92 to ~0.67. The mood is shared by all items, so alpha counts it as consistency; but it is not stress. Alpha is evidence of reliability, not of validity: a scale can be highly consistent about the wrong thing. In real questionnaires, shared noise comes from response styles, such as acquiescence bias, and from items worded so similarly that they measure linguistic overlap.
Exercise 8: PCA of the four scale scores
Run a PCA of the four scale scores (stress, burnout, support, and satisfaction), and interpret the first component from its loadings.
In scikit-learn, pca.components_ contains the principal axes in feature space (transposed loadings). Look at the signs in PC1.
The first component captures 56% of the variation in the four scores. Its loadings are positive for stress (0.55) and burnout (0.53) and negative for support (−0.42) and satisfaction (−0.49): it represents a general dimension of student wellbeing and adjustment, contrasting distress with institutional support and satisfaction. (The overall direction of a component is arbitrary; only the relative contrast between opposite signs matters.)
Exercise 9: Clusters in random data
k-means always returns the number of clusters it is asked for, even when there are none. Run it with four clusters on random uniform data of the same size as the profile data, and compare the share of variation between clusters with the share in the real profiles.
In k-means, inertia_ is the within-cluster sum of squares. Between-cluster SS equals Total SS minus Within-cluster SS.
On purely random uniform data, k-means still produces four clusters capturing ~26% of the total variation: clustering found “groups” where none exist. The real student profiles give ~47%, clearly more, showing that the student data has genuine multidimensional structure. But the comparison also illustrates why clusters must be described as descriptive summaries, not discovered species.
Exercise 10: Why PCA scales the variables
The usarrests.csv dataset contains violent crime and arrest rates across US states in 1973. Run a PCA without scaling and then with scaling, and compare the loadings of the first two components.
Look at the variances first: Assault varies thousands of times more than Murder because of its measurement units. StandardScaler places all features on a standard normal metric.
Without scaling, the first component is almost entirely Assault (loading 1.00) and accounts for ~97% of the total variance simply because assault rates are measured per 100,000 residents and have large numerical variances. After scaling with StandardScaler, every variable counts equally: PC1 (62%) captures overall violent crime (Murder, Assault, Rape), and PC2 (25%) reflects urban population percentage.
Check your understanding
Answer each question in your own words first, then click to see a model answer.
1. Why does a questionnaire use several items for each construct?
Because every item carries the construct plus noise of its own: how the student read the question, their mood, how they use the answer scale. Averaging several items keeps the shared signal and lets the separate noise partly cancel, so the average follows the construct far more closely than any single item.
2. What is the difference between PCA and factor analysis?
PCA summarises the total variation of the variables in a few components, without any theory of where it comes from; it is a data-reduction tool. Factor analysis assumes that latent variables cause the items’ shared variation, and models only that shared part, separating it from each item’s own noise. For checking whether questionnaire items measure the intended constructs, factor analysis is the method that matches the question.
3. An item loads 0.45 on one factor and 0.40 on another. What does this cross-loading mean?
That the item reflects two constructs at once, rather than the one it was written for. Averaged into either scale, it would carry some of the other construct with it. Its wording deserves a look, and it may be better to drop it or to report the results with and without it.
4. A scale has a Cronbach’s alpha of 0.93. Does that show that it is valid?
No. Alpha shows that the items agree with each other, which is reliability. They might all agree about the wrong thing, for example because of a shared response style or very similar wording. A very high alpha can even be a warning sign that the items are near-duplicates. Validity needs other evidence: the factor structure, and correlations with related and unrelated measures.
5. Why should clusters not be described as natural groups of people?
Because clustering always produces the number of groups it is asked for, even in random data, and the result depends on the variables chosen, their scaling, the method, and the number of clusters. Real people differ gradually and overlap. Clusters are useful summaries of typical combinations, and should be reported as such, with the choices that produced them.
Do it yourself
These tasks have no starter code and no answers. Each code space below is empty and ready to run: write your own code, as you would for a thesis. The datasets (questionnaire, items, scores, semesters, profiles, usarrests) and standard packages (numpy, pandas, sklearn, statsmodels) are already imported.
1. Check the burnout scale fully: the correlations between its items, Cronbach’s alpha, and alpha with each item dropped in turn. Decide whether any item weakens the scale, and write the one-sentence report of its reliability.
2. Cluster the students on a different set of variables of your choice (for example, sleep, caffeine, and wellbeing). Choose the number of clusters, describe each cluster by its averages, and argue whether the clusters would be useful to a university.
3. Write the paragraph of a methods section that describes one scale of the questionnaire: the number of items, an example item, the answer scale, how the score was calculated (including any reversed items), and its reliability, with the numbers calculated by code.
4. Run a PCA on usarrests, decide how many components are worth keeping using a scree plot or eigenvalue criterion, and interpret them.
Work on your own computer
The project contains the data and all four parts of this page as a clean Python project: the exercises with blanks, the questions to check your understanding (answers in solutions.py), and the open tasks, each with space to write your code.
- Download chapter09.zip, unzip it, and open the folder in VS Code or your favourite Python editor.
- Install the required packages:
pip install pandas numpy scipy scikit-learn statsmodels matplotlib seaborn