Unsupervised Dimensionality Reduction and Representation

Chapter 23 introduced K-means with a six-patient example, then ran PCA on the full synthetic four-variable table; this chapter starts from that PCA, then adds how to choose the number of components and the nonlinear alternatives. High-dimensional data are difficult to inspect. Dimensionality reduction creates a smaller numerical representation, but “explains variance” is not the same as “predicts disease.” Projection and nonlinear embeddings serve different purposes.

1. PCA: directions of maximum variation

Here is the PCA from chapter 23 again, on the four variables (age, BMI, SBP, pulse). Principal component analysis (PCA) creates new orthogonal axes that capture as much numeric variation as possible, in descending order. It is sensitive to units, outliers, missingness, and input choices.

library(tidymodels)
ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE)
vitals <- ehr |> select(age, bmi, sbp, pulse)
rec <- recipe(~ age + bmi + sbp + pulse, data = vitals) |>
  step_impute_median(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())
scaled <- as.matrix(bake(prep(rec, training = vitals), new_data = vitals))
pca <- prcomp(scaled, center = FALSE, scale. = FALSE)
summary(pca)

Importance of components:
                          PC1    PC2    PC3     PC4
Standard deviation     1.3488 1.0070 0.9929 0.42529
Proportion of Variance 0.4548 0.2535 0.2465 0.04522
Cumulative Proportion  0.4548 0.7083 0.9548 1.00000

round(pca$rotation, 3)   # Loadings: directions of each original variable

         PC1    PC2    PC3    PC4
age   -0.640 -0.372  0.202  0.641
bmi   -0.283  0.876 -0.248  0.303
sbp   -0.705  0.038  0.074 -0.705
pulse -0.117 -0.307 -0.945  0.002

A loading describes how strongly a variable contributes to a component. A score is the location of one observation on a component. The sign of a component can reverse without changing its information content.

plot(pca$x[, 1], pca$x[, 2],
     xlab = "PC1 score", ylab = "PC2 score",
     main = "PCA of synthetic measurements")

↔ Swipe the plot sideways to see all of it.

Do not assume an age/BMI-loaded PC is automatically “diabetes risk.” The outcome was not used to define the PCA axes.

2. Choosing components: scree plot and cumulative variance

The explained variance ratio reports the share of total variance captured by each component. Two common readings use it. The scree plot draws the variance of each component in order; look for the elbow (defined below). The cumulative variance rule keeps enough components to pass a threshold such as 80%. Do not blindly retain “80%” of variance if rare but important directions matter.

share <- pca$sdev^2 / sum(pca$sdev^2)
cbind(component = seq_along(share),
      variance_share = round(share, 3),
      cumulative = round(cumsum(share), 3))

     component variance_share cumulative
[1,]         1          0.455      0.455
[2,]         2          0.254      0.708
[3,]         3          0.246      0.955
[4,]         4          0.045      1.000

Components 1 to 3 hold about 45%, 25% and 25% of the variance; component 4 holds under 5%. Base R draws the scree plot directly from the prcomp object. It shows variances (the squared standard deviations) rather than shares, so the shape is the same.

screeplot(pca, type = "lines", main = "screeplot(pca)")

↔ Swipe the plot sideways to see all of it.

Convention used in this chapter: the elbow is the last component before the line drops steeply. The line falls from component 1 to 2, is nearly flat from 2 to 3, and then drops sharply to component 4, so the elbow suggests 3 components. The 80% cumulative rule also needs 3 (cumulative 0.955, against 0.708 for two).

3. The Kaiser rule

On standardized data every variable has variance 1, so a component that explains less than one original variable adds nothing new. The Kaiser rule keeps components whose eigenvalue is greater than 1. The eigenvalue of a component is its variance, which for prcomp() is sdev^2; the four eigenvalues add up to the number of variables.

eigenvalues <- pca$sdev^2
round(eigenvalues, 2)

[1] 1.82 1.01 0.99 0.18

sum(eigenvalues)

[1] 4

sum(eigenvalues > 1)

[1] 2

The eigenvalues are about 1.82, 1.01, 0.99 and 0.18, and they sum to 4. The rule keeps 2 components, but component 2 (1.01) and component 3 (0.99) are almost the same size and sit on either side of the cut-off. The Kaiser rule is simple but crude and tends to keep too many components when there are many variables; use it as one opinion among several.

4. Parallel analysis

Parallel analysis asks a fairer question: how big would the eigenvalues be if the variables were unrelated noise? Generate random data of the same size (same number of rows and columns), compute its eigenvalues, repeat many times, and average. Keep a real component only if its eigenvalue is larger than the average random eigenvalue at the same position. The rule used here: count components from the first one and stop at the first component that does not beat the random average.

set.seed(26)
n <- nrow(scaled); p <- ncol(scaled)
random_eig <- replicate(200, eigen(cor(matrix(rnorm(n * p), n, p)))$values)
avg_random <- rowMeans(random_eig)
data.frame(component = 1:p,
           real = round(eigenvalues, 2),
           random = round(avg_random, 2),
           keep = eigenvalues > avg_random)

  component real random  keep
1         1 1.82   1.09  TRUE
2         2 1.01   1.02 FALSE
3         3 0.99   0.97  TRUE
4         4 0.18   0.91 FALSE

plot(eigenvalues, type = "b", ylim = c(0, 2),
     xlab = "Component", ylab = "Eigenvalue")
lines(avg_random, type = "b", lty = 2, col = "grey40")
legend("topright", legend = c("real data", "random data (average)"),
       lty = c(1, 2), col = c("black", "grey40"), bty = "n")

↔ Swipe the plot sideways to see all of it.

replicate(200, ...) builds 200 random data sets, and rowMeans() averages their eigenvalues position by position. Only component 1 beats the random line clearly (1.82 against 1.09). Component 2 (1.01 against 1.02) falls just short, so the stop-at-the-first-failure rule gives 1 component. Counting every position separately (the keep column has TRUE, FALSE, TRUE, FALSE) would give 2, because component 3 (0.99 against 0.97) is a near tie. So parallel analysis gives 1, Kaiser gave 2 and the elbow gave 3. These rules often disagree, and that is normal: pick a number you can defend, and say which rule you used. The psych package packages this method as psych::fa.parallel().

5. Choosing by reconstruction error on held-out rows

Another idea is to judge the number of components by how well they rebuild rows the PCA never saw. Fit PCA on training rows only, using the training means and standard deviations. Project the held-out rows onto the first k components, rebuild the four variables from those k scores, and measure the mean squared difference from the true values.

set.seed(26)
raw <- as.matrix(na.omit(vitals))                   # rows with all four values
test_rows <- sample(nrow(raw), 150)
pc_train <- prcomp(raw[-test_rows, ], center = TRUE, scale. = TRUE)   # training means and SDs only
test_z <- scale(raw[test_rows, ], center = pc_train$center, scale = pc_train$scale)

recon_error <- function(k) {
  if (k == 0) return(mean(test_z^2))                    # no components: predict the training mean
  V <- pc_train$rotation[, 1:k, drop = FALSE]
  recon <- test_z %*% V %*% t(V)                        # scores, then back to the 4 variables
  mean((test_z - recon)^2)
}
round(sapply(0:4, recon_error), 3)

[1] 0.941 0.527 0.306 0.041 0.000

The five numbers are the held-out error for 0 to 4 components. This error always falls as k grows, and with all 4 components it is exactly 0, so you cannot pick the minimum. Instead look for where the drop levels off. Here the error falls by about 0.41, 0.22, 0.27 and then only 0.04, so the curve levels off after 3 components, in line with the elbow. That reading rests on the last step alone (0.27 for the third component, then 0.04 for the fourth), so it is a judgment, not a sharp cut. A tool that does pick a number by cross-validation is missMDA::estim_ncpPCA() in the missMDA package; it performs cross-validation.

To use PCA features in a predictive model, put step_pca() inside the recipe fitted on each training fold. Otherwise test-set means, scales, and components leak into training.

6. Reading loadings with a bar chart

Printed loadings are hard to scan. Reshape them to long form with pivot_longer() and draw one bar chart per component with facet_wrap().

loadings <- as.data.frame(pca$rotation)
loadings$variable <- rownames(loadings)
loadings_long <- pivot_longer(loadings, cols = starts_with("PC"),
                              names_to = "component", values_to = "loading")
ggplot(loadings_long, aes(x = loading, y = variable)) +
  geom_col() +
  facet_wrap(~ component, nrow = 1) +
  labs(x = "loading", y = NULL)

↔ Swipe the plot sideways to see all of it.

Each panel is one component and each bar shows how a variable contributes. PC1 is dominated by age and sbp (both negative, so they move together), PC2 by bmi, and PC3 by pulse. Remember that the sign of a whole component is arbitrary: flipping all bars in a panel means the same thing.

7. Biplot: scores and loadings together

A biplot puts the patients (scores) and the variables (loading arrows) on the same PC1-PC2 plane. Points close to an arrow’s direction have high values of that variable, and arrows pointing the same way are correlated variables.

scores <- as.data.frame(pca$x)
arrows <- as.data.frame(pca$rotation[, 1:2])
arrows$variable <- rownames(arrows)
k <- 3                                          # stretch arrows so they are visible
ggplot() +
  geom_point(data = scores, aes(PC1, PC2), colour = "grey60", alpha = 0.5) +
  geom_segment(data = arrows,
               aes(x = 0, y = 0, xend = PC1 * k, yend = PC2 * k),
               arrow = arrow(length = unit(0.2, "cm")), colour = "firebrick") +
  geom_text(data = arrows, aes(PC1 * k * 1.15, PC2 * k * 1.15, label = variable),
            colour = "firebrick")

↔ Swipe the plot sideways to see all of it.

The age and sbp arrows both point to the left, close together, matching their correlation (0.74 for these data; the correlation matrix is printed in chapter 23). bmi points up, mostly along PC2, and the short pulse arrow points down. Base R draws a version with biplot(pca), and the ggfortify package adds autoplot(pca, loadings = TRUE). The arrow length multiplier k is only for display.

8. More than two components

A plane shows two components at a time. To bring a third in, map PC3 to color (or size) on the PC1-PC2 scatter.

ggplot(scores, aes(PC1, PC2, colour = PC3)) +
  geom_point(alpha = 0.7) +
  scale_colour_viridis_c()

↔ Swipe the plot sideways to see all of it.

In this plot the PC3 colors are mixed all over the PC1-PC2 plane, with no visible gradient in any direction, so PC3 adds no visible pattern to this view. Interactive 3D scatter plots exist, for example through the plotly package, and let you rotate the cloud, but a static plot with color is easier to put in a report.

9. SVD and matrix factorization

Singular value decomposition (SVD) factorizes a matrix into orthogonal directions and strength weights. It underlies PCA and many recommendation/text methods. For a centered matrix, the leading SVD directions relate to PCA directions.

s <- svd(scaled)
round(s$d[1:3], 3)

[1] 33.011 24.646 24.301

dim(s$u)

[1] 600   4

dim(s$v)

[1] 4 4

Example: the singular values are on a different scale from PCA’s standard deviations, but they carry the same information. Divide them by sqrt(n - 1) and they match pca$sdev.

round(s$d / sqrt(nrow(scaled) - 1), 3)

[1] 1.349 1.007 0.993 0.425

round(pca$sdev, 3)

[1] 1.349 1.007 0.993 0.425

The two lines are identical (1.349, 1.007, 0.993, 0.425). s$v holds the loadings (4 rows, one per variable) and s$u the scaled scores (600 rows, one per patient), which is why dim(s$u) printed 600 rows.

NMF (nonnegative matrix factorization) decomposes nonnegative data into additive parts; useful for counts, topic components, and nonnegative signal patterns. R: NMF package. It assumes data and factors have the required nonnegative structure.

10. Factor analysis and ICA

Factor analysis models covariance as arising from latent variables plus measurement-specific noise; it is not the same objective as PCA. R: stats::factanal().

Example: fit one latent factor to the four standardized measurements.

fa <- factanal(scaled, factors = 1)
round(fa$uniquenesses, 3)

  age   bmi   sbp pulse 
0.451 0.886 0.005 0.995 

signif(fa$PVAL, 2)                    # p-value of the test that one factor is enough

objective 
  5.1e-26 

print(fa$loadings, cutoff = 0.3)

Loadings:
      Factor1
age   0.741  
bmi   0.338  
sbp   0.998  
pulse        

               Factor1
SS loadings      1.663
Proportion Var   0.416

The loading of sbp on the factor is about 1 and its uniqueness (the part not shared with the factor) is near 0, so the factor is mostly systolic pressure; age loads at 0.74 and bmi weakly. pulse has uniqueness near 1, so it shares almost nothing with the others and the printout hides its loading under the 0.3 cutoff. The model is a one-factor story for four variables, and the printed p-value (fa$PVAL, about 5e-26) is far below 0.05, so the test rejects one factor as sufficient. Treat it as a demonstration of the output, not a finding.

Independent component analysis (ICA) seeks statistically independent latent signals, often used in signal processing. R: fastICA. Independence assumptions and component identifiability require care.

Example scenario for ICA: two people talk at once and two microphones each record a mix of both voices. PCA would return directions that explain the most variance, which are still mixtures. ICA looks for two output signals that are statistically independent of each other, and for this problem that recovers one voice per signal. Health data rarely fit such a clean story, which is why ICA appears mostly in EEG and other signal work.

11. Nonlinear geometry: t-SNE and UMAP

PCA is linear. t-SNE and UMAP construct nonlinear embeddings, most commonly for exploratory visualization of complex high-dimensional data.

t-SNE/UMAP plots depend on parameters, seeds, densities, and preprocessing. Do not compare the physical distance between two plotted clusters as though it were a clinical risk difference. They also require clear out-of-sample handling for predictive use.

Example scenario: a t-SNE plot of 600 patients shows three islands. Re-running it with a different seed and perplexity (the neighborhood size) gives two islands and a thin bridge between them. Both plots come from the same data, so the number of islands is a property of the settings as much as of the patients. Before writing “three subtypes”, check the islands against something outside the plot, such as the original measurements, and repeat the embedding with several seeds. A PCA of the same data gives the same axes every time, which is why it is the safer first look.

Optional demonstration (additional packages, not required for core CI):

# install.packages(c("Rtsne", "uwot"))
# ts <- Rtsne::Rtsne(scaled, perplexity = 20, check_duplicates = FALSE)
# um <- uwot::umap(scaled, n_neighbors = 15, n_components = 2)

12. Sparse and unstructured data

Text documents, diagnosis-code count tables, images, and sensor time series need different representations. Frequent approaches include SVD/LSA for sparse text matrices, topic models, embeddings from language models, CNN feature extractors, and sequence embeddings. Before extracting features from clinical notes, obtain approved privacy/security review; never paste real notes into an unapproved service.

Example: a table of diagnosis-code counts per patient is mostly zeros. Here are six made-up patients and four codes. SVD finds the groups of codes that appear together.

codes <- matrix(c(2,0,1,0,  3,0,2,0,  0,2,0,3,  0,1,0,2,  1,0,1,0,  0,3,0,1),
                nrow = 6, byrow = TRUE,
                dimnames = list(paste0("patient", 1:6),
                                c("hypertension", "statin", "diabetes", "insulin")))
lsa <- svd(codes)
round(lsa$d, 2)

[1] 5.00 4.46 1.73 0.39

round(lsa$v[, 1:2], 2)

      [,1]  [,2]
[1,]  0.00 -0.84
[2,] -0.71  0.00
[3,]  0.00 -0.54
[4,] -0.71  0.00

The first two singular values (5.00 and 4.46) are much larger than the others, so two directions describe most of the table. Their loadings in lsa$v show what they are: the first pairs statin with insulin, and the second pairs hypertension with diabetes. The signs of the printed columns are arbitrary. Real text or code tables have thousands of columns, and the same calculation (run on a sparse matrix) is called LSA.

13. Validation for representation methods

For description, examine reconstruction/variance, robustness across seeds, and whether features reflect device/site artifacts. For supervised prediction, learn representation transformations using only training data, evaluate downstream outcomes in appropriate held-out groups, and compare to the original features. Prevent identifiers and post-outcome signals from entering representation learning.

Example: learn the PCA on training rows only, then apply it unchanged to held-out rows. In a recipe this is what prep() on the training data and bake() on new data do.

set.seed(26)
train_rows <- sample(nrow(vitals), 450)
rec_pca <- recipe(~ age + bmi + sbp + pulse, data = vitals[train_rows, ]) |>
  step_impute_median(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors()) |>
  step_pca(all_numeric_predictors(), num_comp = 2)
prepped <- prep(rec_pca)                       # medians, means, SDs and loadings come from 450 training rows
bake(prepped, new_data = vitals[-train_rows, ]) |> head(3)

# A tibble: 3 x 2
     PC1     PC2
   <dbl>   <dbl>
1 -1.80  -0.547 
2  0.517  0.0198
3  1.29   0.206 

The three printed rows are held-out patients expressed on PC1 and PC2, using only what the 450 training rows taught the recipe. Fitting the PCA on all 600 rows first would let the held-out rows shape the axes, which is the leak this section warns about.

Try it

1. How many components are required for at least 80% variance in the toy vitals?

Show solution to question 1

Answer: 3 components. Two components hold 70.8% of the variance and three hold 95.5%, so 3 is the smallest number that reaches 80%.

Step by step

1. The chapter already computed share, the variance share of each component. Take its running total with cumsum().
2. Find the first position where the total reaches 0.80.

round(cumsum(share), 3)

Output

[1] 0.455 0.708 0.955 1.000

which(cumsum(share) >= 0.80)[1]

Output

[1] 3

3. The running totals are 0.455, 0.708, 0.955 and 1. The first one at or above 0.80 is the third, so which() returns 3. This depends on the exact columns and on median imputation, and it says nothing about whether the components predict anything.

2. Remove pulse, repeat PCA, and compare loadings and variance shares.

Show solution to question 2

Answer: Without pulse the three components hold 60.3%, 33.7% and 6.0% of the variance (with pulse they were 45.5%, 25.4%, 24.6% and 4.5%). PC1 is still mostly age and sbp, and PC2 is still bmi; only the sign of PC2 flipped, which carries no meaning.

Step by step

1. Take the three standardized columns from the chapter’s scaled matrix and run prcomp() on them.
2. Print the loadings and the variance shares.

scaled3 <- scaled[, c("age", "bmi", "sbp")]
pca3 <- prcomp(scaled3, center = FALSE, scale. = FALSE)
round(pca3$rotation, 3)

Output

       PC1    PC2    PC3
age -0.644  0.417  0.642
bmi -0.286 -0.909  0.303
sbp -0.710 -0.012 -0.705

round(pca3$sdev^2 / sum(pca3$sdev^2), 3)

Output

[1] 0.603 0.337 0.060

round(pca$rotation[, 1:3], 3)       # the four-variable loadings, for comparison

Output

         PC1    PC2    PC3
age   -0.640 -0.372  0.202
bmi   -0.283  0.876 -0.248
sbp   -0.705  0.038  0.074
pulse -0.117 -0.307 -0.945

3. Compare. PC1 is -0.644 for age and -0.710 for sbp here, against -0.640 and -0.705 before. PC2 is bmi in both (-0.909 here, +0.876 before), so the sign flipped but the meaning did not. The four-variable PC3 was almost all pulse; that component is gone, which is why the first two shares grow: the total variance is now 3 instead of 4.

3. State one fundamental difference between PCA and t-SNE.

Show solution to question 3

Answer: PCA is a linear projection with fixed, interpretable axes (the loadings) and the same result on every run. t-SNE is a nonlinear embedding that keeps local neighborhoods, has no loadings, depends on the seed and perplexity, and distorts global distances.

Step by step

1. PCA builds each new axis as a weighted sum of the original variables, so you can read what an axis means and project new rows onto it.
2. t-SNE only arranges points in 2D so that near neighbors stay near. The plot coordinates have no meaning of their own.
3. Because of that, a gap between two t-SNE islands is not evidence of two diseases, and the plot is mainly for exploration.

4. Explain the difference between unsupervised representation learning and supervised feature selection.

Show solution to question 4

Answer: Unsupervised representation learning (PCA, for example) builds new features from the predictors alone, without looking at the outcome. Supervised feature selection uses the outcome to decide which features to keep.

Step by step

1. PCA used only age, bmi, sbp and pulse; diabetes never entered, which is why an age/BMI axis is not a “diabetes risk” axis.
2. Supervised selection (keeping the variables most related to diabetes) uses the labels, so a choice made on all rows can overfit to them.
3. Each step must therefore be learned inside the training fold of a validation scheme, and the supervised one needs this more urgently because it has seen the outcome.

5. Why is a full-data PCA inappropriate if it will later feed a locked-test predictor?

Show solution to question 5

Answer: The full-data PCA computes its means, scales and loadings partly from the test rows, so information from the locked test set leaks into the training features and the test result can look better than it should be (for an unsupervised method like PCA the leak is usually small, but it is still a leak).

Step by step

1. A locked test set is meant to behave like future patients nobody has seen.
2. A PCA fitted on all rows lets those rows shape the axes, so the model is trained on features that already contain test-set information.
3. The fix is to fit the recipe (impute, normalize, step_pca()) on the training rows, or on each training fold during cross-validation, and only apply it to the test rows, as the prep() / bake() example in section 13 does.

6. Apply the Kaiser rule and parallel analysis to the toy vitals with pulse removed. Do they agree?

Show solution to question 6

Answer: Yes. Kaiser keeps 2 components (eigenvalues 1.81 and 1.01 are above 1) and parallel analysis, counted until the first failure, also keeps 2. Component 2 (1.01 against a random average of 1.00) only just passes, so treat 2 as a close call.

Step by step

1. Take the eigenvalues of the three-variable PCA from item 2 and count those above 1 (Kaiser).
2. Simulate random data of the same size, average its eigenvalues and compare position by position, stopping at the first failure.

ev3 <- pca3$sdev^2
set.seed(26)
n3 <- nrow(scaled3)
random3 <- rowMeans(replicate(200, eigen(cor(matrix(rnorm(n3 * 3), n3, 3)))$values))
list(eigenvalues = round(ev3, 2),
     random = round(random3, 2),
     kaiser = sum(ev3 > 1),
     parallel = which(!(ev3 > random3))[1] - 1)

Output

$eigenvalues
[1] 1.81 1.01 0.18

$random
[1] 1.06 1.00 0.94

$kaiser
[1] 2

$parallel
[1] 2

3. The eigenvalues are 1.81, 1.01 and 0.18 against random averages of 1.06, 1.00 and 0.94. Kaiser counts 2. In parallel analysis the first component that fails to beat the random average is the third, so the count is 2 as well. The margin on component 2 is about 0.01, so a different seed could move it. With pulse included (section 4) the same stop-at-the-first-failure rule kept only 1 component, because component 2 fell just short (1.01 against 1.02). Here it just passes (1.01 against 1.00), so the count depends on which variables go in, and component 2 is a near tie in both versions.

7. Draw the loadings bar chart for the PCA without pulse.

Show solution to question 7

Answer: PC1 is carried by age and sbp (-0.644 and -0.710), PC2 by bmi (-0.909), and PC3 contrasts age (0.642) with sbp (-0.705).

Step by step

1. Turn the loadings of pca3 into a data frame and reshape it to long form with pivot_longer().
2. Draw one panel per component with facet_wrap().

load3 <- as.data.frame(pca3$rotation)
load3$variable <- rownames(load3)
ggplot(pivot_longer(load3, starts_with("PC"), names_to = "component", values_to = "loading"),
       aes(x = loading, y = variable)) +
  geom_col() +
  facet_wrap(~ component, nrow = 1) +
  labs(x = "loading", y = NULL)

↔ Swipe the plot sideways to see all of it.

3. Read the panels against the printed table in item 2. The long bars in each panel are the variables that define that component. The sign of a whole panel is arbitrary.

8. Make a biplot of PC1 against PC3 instead of PC1 against PC2, and say which variable dominates PC3.

Show solution to question 8

Answer: pulse dominates PC3: its loading is -0.945, against 0.202 for age, -0.248 for bmi and 0.074 for sbp. In the plot the pulse arrow is long and points along the PC3 axis.

Step by step

1. Print the PC3 loadings of the four-variable PCA.
2. Build the arrows from the PC1 and PC3 columns of the loadings and plot them over the PC1 and PC3 scores.

round(pca$rotation[, "PC3"], 3)

Output

   age    bmi    sbp  pulse 
 0.202 -0.248  0.074 -0.945 

arrows13 <- as.data.frame(pca$rotation[, c("PC1", "PC3")])
arrows13$variable <- rownames(arrows13)
ggplot() +
  geom_point(data = scores, aes(PC1, PC3), colour = "grey60", alpha = 0.5) +
  geom_segment(data = arrows13,
               aes(x = 0, y = 0, xend = PC1 * 3, yend = PC3 * 3),
               arrow = arrow(length = unit(0.2, "cm")), colour = "firebrick") +
  geom_text(data = arrows13, aes(PC1 * 3.4, PC3 * 3.4, label = variable),
            colour = "firebrick")

↔ Swipe the plot sideways to see all of it.

3. The printed loadings show one large value (-0.945 for pulse), so PC3 is essentially a pulse axis. This agrees with the bar chart in section 6, and with the PC3 colors in section 8 showing no pattern on the PC1-PC2 plane.

9. Re-run the held-out reconstruction error with a different seed. Does the point where the drop levels off move?

Show solution to question 9

Answer: No. With seed 1 the errors for 0 to 4 components are 0.929, 0.497, 0.272, 0.039 and 0, and the drops are 0.432, 0.225, 0.233 and 0.039, so the curve still levels off after 3 components. The values differ slightly from seed 26 (0.941, 0.527, 0.306, 0.041, 0).

Step by step

1. Wrap the chapter’s steps in a function of the seed: pick 150 test rows, fit PCA on the rest, project the test rows and measure the error for 0 to 4 components.
2. Run it for seed 26 and seed 1 and print the differences between neighbours.

recon_seed <- function(seed) {
  set.seed(seed)
  rows <- sample(nrow(raw), 150)
  pt <- prcomp(raw[-rows, ], center = TRUE, scale. = TRUE)
  z <- scale(raw[rows, ], center = pt$center, scale = pt$scale)
  sapply(0:4, function(k) {
    if (k == 0) return(mean(z^2))
    V <- pt$rotation[, 1:k, drop = FALSE]
    mean((z - z %*% V %*% t(V))^2)
  })
}
errs <- rbind(seed_26 = recon_seed(26), seed_1 = recon_seed(1))
colnames(errs) <- paste0("k", 0:4)
round(errs, 3)

Output

           k0    k1    k2    k3 k4
seed_26 0.941 0.527 0.306 0.041  0
seed_1  0.929 0.497 0.272 0.039  0

round(-t(apply(errs, 1, diff)), 3)      # how much each extra component lowers the error

Output

           k1    k2    k3    k4
seed_26 0.414 0.220 0.265 0.041
seed_1  0.432 0.225 0.233 0.039

3. Both rows fall with each added component and end at 0. In both, the first three components each lower the error by 0.2 to 0.4 and the fourth by only about 0.04, so the leveling off is after 3. The exact numbers move with the split, which is why this is a judgment and not a sharp cut.

10. Color the PC1-PC2 scatter by the original age instead of PC3. Why does it show a clear pattern when PC3 did not?

Show solution to question 10

Answer: The colors form a clear gradient along PC1, because age has a large loading on PC1 (-0.64) and the correlation between standardized age and the PC1 score is -0.86. PC3 is mostly pulse and is only weakly related to the other variables (its correlation with age is 0.20), so it showed no pattern on the PC1-PC2 plane.

Step by step

1. Add the standardized age column to the scores and map it to color.

ggplot(cbind(scores, age = scaled[, "age"]), aes(PC1, PC2, colour = age)) +
  geom_point(alpha = 0.7) +
  scale_colour_viridis_c()

↔ Swipe the plot sideways to see all of it.

2. Check the link with numbers.

round(c(with_PC1 = cor(scaled[, "age"], scores$PC1),
        with_PC3 = cor(scaled[, "age"], scores$PC3)), 2)

Output

with_PC1 with_PC3 
   -0.86     0.20 

3. Age is strongly tied to PC1 (-0.86, so older patients have lower PC1 scores) and barely to PC3 (0.20). Color shows a pattern when the colored variable lines up with an axis of the plot, and nothing when it does not.