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 some of R’s built-in 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 functions that come with R: prcomp() for PCA, factanal() for factor analysis, kmeans() and hclust() for clustering, and a few lines of code for Cronbach’s alpha, which you build yourself in Exercise 3. The chapter’s psych and factoextra packages give the same results with more detail on your own computer. 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.
Base R’s function for principal component analysis is named after it, with “r” in front for R.
support_pca <- prcomp(support_items, scale. = TRUE)
summary(support_pca)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. (na.omit() 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.
factors is the number of factors to extract. rotation = "promax" lets the factors correlate, like the chapter’s oblimin rotation.
fa3 <- factanal(na.omit(items), factors = 3, rotation = "promax")
print(fa3$loadings, cutoff = 0.3, sort = TRUE)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 factanal() uses only the 396 students who answered every item, while the chapter’s fa() used everyone.
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) \]
apply(x, 2, f) applies the function f to every column (the 2) of x; here f is the function for the variance. On a 1-to-5 scale, subtracting an answer from 6 reverses it.
cronbach <- function(x) {
x <- na.omit(x)
k <- ncol(x)
k / (k - 1) * (1 - sum(apply(x, 2, var)) / var(rowSums(x)))
}
stress_reversed <- stress_items |> mutate(stress_4 = 6 - stress_4)Without reversing, alpha is 0.57, below the usual threshold of 0.7; 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. (psych’s alpha() gives the same values, and warns that stress_4 probably should be reversed.)
Exercise 4: Three clusters
Run k-means with 3 clusters on the profile data, describe the clusters, and identify which profiles from the 4-cluster solution were merged.
centers is the number of clusters. In the final table, a column of the four-cluster solution whose students all fall into one row has been absorbed into that three-cluster group.
set.seed(123)
k3 <- kmeans(profile_data, centers = 3, nstart = 25)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 table 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. (Cluster numbers are arbitrary and can differ from run to run; read the averages, not the numbers.)
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 a clear 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 research [reference], and 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 groups.” The paragraph names the criteria, admits that they were not decisive, gives the substantive reason for the choice, and says how the clusters should be read.
Exercise 6: Noisier items
Change the simulation at the start of the chapter so that each item has more noise (sd = 2), 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 standard deviation of rnorm(); the chapter used 1.
simulated_items <- sapply(1:20, function(i) true_stress + rnorm(600, sd = 2))With twice the noise, one item correlates only 0.41 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 a vector of 600 values to the 600 × 6 matrix adds each student’s mood to all six of their items.
items_mood <- items_clean + moodWith 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 a tendency to agree with every statement, and from items worded so similarly that they measure the wording.
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.
scale. = TRUE standardises each variable before the PCA, so that each counts equally. The loadings are in the rotation part of the result: look at their signs in the first column.
scores_pca <- prcomp(scale_scores, scale. = TRUE)
round(scores_pca$rotation, 2)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 is a general dimension of how a student is doing, with distress at one end and support and satisfaction at the other. (The overall sign of a component is arbitrary; only the contrast between the two pairs matters.) The four scales are distinct, as the factor analysis showed, but they share a large common core.
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 data of the same size as the profile data, and compare the share of variation between clusters with the share in the real profiles.
betweenss is the variation between the cluster centres; divided by the total sum of squares, it gives the share of all variation that the clusters capture. The name of the total follows the same pattern as the real line below it.
c(random = k_random$betweenss / k_random$totss,
real = k_real$betweenss / k_real$totss)On purely random data, k-means still produced four tidy clusters of 146 to 157 points each, and they capture 26% of the variation: clustering found “groups” where none exist. The real profiles give 47%, clearly more, so the student data has more structure than random noise. But the comparison also shows why clusters must be described as summaries, not discovered types: a large part of what any clustering captures would appear even in data without groups.
Exercise 10: Why PCA scales the variables
R’s built-in USArrests data gives four crime-related rates for the 50 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: one variable varies thousands of times more than another, only because of its units.
scaled <- prcomp(USArrests, scale. = TRUE)Without scaling, the first component is simply Assault (a loading of 1.00) and “explains” 97% of the variation, because assault rates are counted per 100,000 people and vary more than 300 times as much as murder rates. After scaling, every variable counts equally. The first component (62%) loads on murder, assault, and rape together: an overall level of violent crime. The second (25%) is mainly the urban population. Variables measured in different units should almost always be scaled before a PCA; items on the same answer scale, as in a questionnaire, need it less.
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 data (questionnaire, items, scores, semesters, profiles, and profile_data) and the dplyr package are already loaded, and R’s built-in datasets are always available.
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, first-semester 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. R’s built-in attitude data records seven ratings of their supervisors by the employees of 30 departments of a large company. Run a PCA, decide how many components are worth keeping, and interpret them.
Work on your own computer
The project contains the data and all four parts of this page as an R script: the exercises with blanks, the questions to check your understanding (answers in solutions.R), and the open tasks, each with space to write your code. Its solutions also show the chapter’s psych functions, fa() and alpha(), side by side with the base R versions.
- Download chapter09.zip, unzip it, and double-click
chapter09.Rproj. - Or type this one line in RStudio’s Console:
usethis::use_course("https://polla-fattah.github.io/data2thesis_r/playground/chapter09.zip")