library(tidyverse)
library(vegan)
library(ggrepel)
library(fastDummies)Activity 2.2 - multidimensional scaling
Plot hint
Make sure your plots are big enough by specifying:
#| fig-height: 9
#| fig-width: 9in your R chunks.
SUBMISSION INSTRUCTIONS
Submit both:
- Your .qmd;
- A link to your rendered html.
Relevant packages (not exhaustive!)
Question 1
Suppose we have 3 data points in 2 dimensions:
A)
Using taxicab as the distance metric, find all 3 pairwise distances between points in the full dimension space.
A-B = 4 A-C = 6 B-C = 4
B)
Consider three one-dimensional ordinations of this space. In the first ordination, the points are collapsed to their vertical coordinates; in the second, they are collapsed to their horizontal coordinates; and in the third, they are collapsed vertically to the y=x line:
Find the stress of each of these ordinations. Which best represents the full-dimensional distances between points?
AB = sqrt(-1^2 + 3^2) = 3.162 AC = sqrt(2^2 + 4^2) = 4.472 BC = sqrt(3^2 + 1^2) = 3.162
1 = ((3.162 - 3)^2 + (4.472 - 4)^2 + (3.162 - 1)^2) = 2.219 2 = ((3.162 - 1)^2 + (4.472 - 2)^2 + (3.162 - 3)^2) = 3.288 3 = ((3.162 - 1.414)^2 + (4.472 - 2.828)^2 + (3.162 - 4.243)^2) = 2.632
1 is the best representation
Question 2
Consider the following MDS biplot of a 2-dimensional representation of 4-dimensional data:
Which of these most likely represents the correlation matrix between the four variables?
## Matrix 1:
A B C D
A 1.00
B 0.90 1.00
C 1.00 0.39 1.00
D 0.39 -0.06 0.06 1.00## Matrix 2:
A B C D
A 1.00
B -1.00 1.00
C 0.94 -0.94 1.00
D 0.06 0.39 -0.06 1.00## Matrix 3:
A B C D
A 1.00
B -1.00 1.00
C -0.94 0.94 1.00
D 0.06 -0.06 0.39 1.00## Matrix 4:
A B C D
A 1.00
B -1.00 1.00
C 0.94 -0.94 1.00
D 0.06 -0.06 0.39 1.00MATRIX 4
Question 3
Consider the correlation matrix between three variables given by:
A B C
A 1.00 0.08 -0.14
B 0.08 1.00 0.98
C -0.14 0.98 1.00Now consider four biplots below showing the 3-dimensional data set reduced to 2 dimensions using NMDS:
Which plot does the correlation matrix belong to?
PLOT 4
Question 4
The data for this problem come from Kaggle. We have various numeric characterstics of 777 universities, both public and private. See the Kaggle page for the variable descriptions.
college <- read.csv('Data/College.csv')The goal is to understand how universities differ with respect to the numeric variables measured, and how these characteristics compare between private and public universities.
A)
What is the best way to measure dissimilarities between universities with respect to the measured variables? Is scaling important before measuring dissimilarity? Why/why not?
Use euclidian distance since all the variables are quantitative. scaling is important because the variables have different ranges, so bigger values like Apps would have more influence on the distance. standardizing makes sure all variables contribute more equally.
B)
Perform non-metric multidimensional scaling using a 2-dimensional ordination, using the dissimilarity measure you chose in A). Add the ordination coordinates to the college data frame. What is the stress of the solution, and what does this indicate? Explain as if to someone who knows nothing about NMDS.
college_scaled <- scale(college[sapply(college, is.numeric)])
nmds <- metaMDS(college_scaled,
distance = "euclidean",
k = 2,
autotransform = FALSE)'comm' has negative data: 'autotransform', 'noshare' and 'wascores' set to FALSERun 0 stress 0.1295778
Run 1 stress 0.1298827
... Procrustes: rmse 0.005785367 max resid 0.1458325
Run 2 stress 0.1297293
... Procrustes: rmse 0.007976541 max resid 0.1542352
Run 3 stress 0.1313113
Run 4 stress 0.132773
Run 5 stress 0.1422159
Run 6 stress 0.1292592
... New best solution
... Procrustes: rmse 0.007328955 max resid 0.1546231
Run 7 stress 0.1291599
... New best solution
... Procrustes: rmse 0.00698521 max resid 0.1836877
Run 8 stress 0.1552981
Run 9 stress 0.1296258
... Procrustes: rmse 0.00697557 max resid 0.1836723
Run 10 stress 0.1313787
Run 11 stress 0.1362616
Run 12 stress 0.1310413
Run 13 stress 0.13033
Run 14 stress 0.1290051
... New best solution
... Procrustes: rmse 0.004772747 max resid 0.1184581
Run 15 stress 0.1302863
Run 16 stress 0.1310625
Run 17 stress 0.1311453
Run 18 stress 0.1293495
... Procrustes: rmse 0.004544116 max resid 0.1181788
Run 19 stress 0.1376565
Run 20 stress 0.1299807
*** Best solution was not repeated -- monoMDS stopping criteria:
3: no. of iterations >= maxit
3: stress ratio > sratmax
14: scale factor of the gradient < sfgrmincollege$NMDS1 <- scores(nmds)[,1]
college$NMDS2 <- scores(nmds)[,2]
nmds$stress[1] 0.1290051The stress is 0.1293, which means the 2D plot does a pretty good job of showing the differences between colleges. There is some distortion, but the plot is still a reasonable representation of the distances between colleges.
C)
Create a stress plot “from scratch,” by creating a data frame consisting of 2 columns:
Column 1: a column of the pairwise distances between the colleges in full-dimensional space;
Column 2: a column of the pairwise Euclidean distances between the colleges in the new 2-dimensional space.
(Hint: applying as.vector() to a vegdist object will form a vector of the \(\frac{n(n-1)}{2}\) unique distances.)
Then use ggplot to plot these two distances versus each other (hint: use shape='.' to speed up the plot rendering - there are thousands of pairwise distances to plot!) Add the value of the correlation between these distances to your plot title. Is the size of this correlation a good thing or a bad thing?
full_dist <- as.vector(vegdist(college_scaled, method = "euclidean"))
twoD_dist <- as.vector(dist(scores(nmds)))
stress_df <- data.frame(full_dist, twoD_dist)
r <- cor(full_dist, twoD_dist)
ggplot(stress_df, aes(x = full_dist, y = twoD_dist)) +
geom_point(shape = ".") +
labs(
title = paste("Correlation =", round(r, 3), ")"),
x = "Distances in Full-Dimensional Space",
y = "Distances in 2D Space"
) +
theme_minimal()
A high positive correlation is good because it means the distances in the 2D plot are similar to the original distances between colleges. A low correlation means the 2D plot does a worse job representing those distances. However, NMDS focuses more on preserving the rank order of distances than matching their exact values.
D)
Use envfit() to fit all numeric variables to the MDS coordinates. Extract the arrow vector coordinates and \(R^2\) and put them in their own data frame. Sort this data set by \(R^2\). Which variables are most and least responsible for separating points in this lower dimensional space?
fit <- envfit(nmds, college_scaled, permutations = 999)
arrows <- scores(fit, display = "vectors")
envfit_df <- data.frame(
Variable = rownames(arrows),
NMDS1 = arrows[, 1],
NMDS2 = arrows[, 2],
R2 = fit$vectors$r
) %>%
arrange(desc(R2))
envfit_df Variable NMDS1 NMDS2 R2
F.Undergrad F.Undergrad 0.778565109 -0.35729378 0.733822478
Apps Apps 0.852914837 -0.04135915 0.729174299
Accept Accept 0.833663441 -0.18397478 0.728841451
Enroll Enroll 0.800010688 -0.29334845 0.726070415
Top10perc Top10perc 0.414714302 0.71085405 0.677301435
Outstate Outstate 0.099306062 0.80506924 0.657998177
Expend Expend 0.274592049 0.74426497 0.629331136
Top25perc Top25perc 0.456344856 0.62746829 0.601967088
PhD PhD 0.633225892 0.42133391 0.578497297
Terminal Terminal 0.579151874 0.43418339 0.523932111
P.Undergrad P.Undergrad 0.482073141 -0.53762941 0.521439898
Grad.Rate Grad.Rate 0.139212658 0.65700671 0.451037977
perc.alumni perc.alumni 0.005708432 0.64811260 0.420082533
S.F.Ratio S.F.Ratio 0.043772070 -0.63809237 0.409077868
Room.Board Room.Board 0.161890811 0.58844160 0.372472149
Personal Personal 0.224416155 -0.38887633 0.201587412
Books Books 0.031927476 0.06375265 0.005083765Most responsable: Apps, F. Undergrad Least responsable: Books, Personal
E)
Use ggplot to plot the 2 dimensional ordination and the labeled arrows. Color-code the points by public/private; optionally, add 2SD ellipses for each college type as well. Which variable(s) are primarily separating samples along the horizontal axis? Vertical axis? What does this plot tell you about how private/public universities differ with respect to the measured variables?
mult <- 0.8 * max(abs(c(college$NMDS1, college$NMDS2))) /
max(sqrt(envfit_df$NMDS1^2 + envfit_df$NMDS2^2))
envfit_df <- envfit_df %>%
mutate(xend = NMDS1 * mult, yend = NMDS2 * mult)
ggplot(college, aes(x = NMDS1, y = NMDS2)) +
geom_point(aes(color = Private), alpha = 0.5, size = 1.5) +
geom_segment(data = envfit_df,
aes(x = 0, y = 0, xend = xend, yend = yend),
arrow = arrow(length = unit(0.2, "cm")),
inherit.aes = FALSE) +
geom_text(data = envfit_df,
aes(x = xend * 1.1, y = yend * 1.1, label = Variable),
size = 3, inherit.aes = FALSE) +
labs(title = "NMDS of US Colleges",
color = "Private") +
coord_equal() +
theme_classic()
F)
For each of the most important variables, create side-by-side boxplots that you can use to compare these variables for public and private universities. Comment on how these plots corroborate your NMDS biplot from the previous question.
top_vars <- envfit_df %>%
slice_max(R2, n = 6) %>%
pull(Variable)
college %>%
select(Private, all_of(top_vars)) %>%
pivot_longer(-Private, names_to = "Variable", values_to = "Value") %>%
ggplot(aes(x = Private, y = Value, fill = Private)) +
geom_boxplot() +
facet_wrap(~ Variable, scales = "free_y") +
labs(title = "Most Important Variables: Public vs. Private",
x = "Private", y = "Value") +
theme_minimal() +
theme(legend.position = "none")
Question 5
For this problem we will analyze the mushroom data, which are various categorical characteristics of poisonous and edible mushrooms. The variables used to measure distance are all except the class variable. See the data set Kaggle site for variable descriptions. Note that stalk.root has several unknown values indicated by ?, which we will remove. We’ll also de-select some variables that are mostly constant and sample to only 1000 mushrooms for simplicity.
set.seed(1112)
mushroomsc <- (read.csv('Data/mushrooms-with-columnnames.csv')
%>% filter(stalk.root!="?")
%>% dplyr::select(-veil.type, -veil.color,-gill.attachment)
%>% sample_n(1000)
)A)
Create an appropriate dissimilarity matrix for this data set, using all variables except class.
Gower
library(cluster)
mushroom_vars <- mushroomsc %>%
dplyr::select(-class) %>%
mutate(across(everything(), as.factor))
mushroom_dist <- daisy(mushroom_vars, metric = "gower")
as.matrix(mushroom_dist)[1:5, 1:5] 1 2 3 4 5
1 0.0000000 0.4736842 0.5789474 0.6315789 0.2105263
2 0.4736842 0.0000000 0.1578947 0.5789474 0.4736842
3 0.5789474 0.1578947 0.0000000 0.6315789 0.4736842
4 0.6315789 0.5789474 0.6315789 0.0000000 0.5789474
5 0.2105263 0.4736842 0.4736842 0.5789474 0.0000000B)
Use the dissimilarity matrix to perform NMDS. (NOTE: by default, metaMDS algorithm starts from 20 different random starts to try to make sure the optimal solution is found. This is a good idea in practice, but with large n, even of 1000, each try takes a while, so 20 will take forever! For now, set trymax=1 to speed things up.) Is there much difference in stress between a 2 and 3 dimensional configuration?
set.seed(1112)
nmds2 <- metaMDS(as.dist(mushroom_dist), k = 2, trymax = 1)Run 0 stress 0.1439484
Run 1 stress 0.1476318
*** Best solution was not repeated -- monoMDS stopping criteria:
1: scale factor of the gradient < sfgrminset.seed(1112)
nmds3 <- metaMDS(as.dist(mushroom_dist), k = 3, trymax = 1)Run 0 stress 0.09984306
Run 1 stress 0.09974373
... New best solution
... Procrustes: rmse 0.002823357 max resid 0.05574656
*** Best solution was not repeated -- monoMDS stopping criteria:
1: scale factor of the gradient < sfgrminc(stress_2D = nmds2$stress, stress_3D = nmds3$stress) stress_2D stress_3D
0.14394843 0.09974373 The 2D stress is 0.144 and the 3D stress is 0.100, so adding another dimension makes the fit better. The 2D plot is still decent and easier to understand, but it does miss some of the differences between mushrooms.
C)
Focus on the 2-dimensional ordination. Use envfit to obtain relationships between the mushroom characteristics and these MDS coordinates. (NOTE: since the variables are factors here, rather than regressing the variables on the coordinates envfit() simply averages the coordinates for each level of the factors; these are contained in the $factors$centroids object.) Create a data frame of these centroids. Then create a biplot with the following specifications:
- The MDS coordinates on horizontal and vertical;
- Points color-coded by poisonous/edibility;
- Points of different shape and color for the centroids;
- Labels for the centroids with measures taken to reduce overlapping labels
Is this plot readable?
mushroomsc$NMDS1 <- scores(nmds2)[, 1]
mushroomsc$NMDS2 <- scores(nmds2)[, 2]
fit_m <- envfit(nmds2, mushroom_vars, permutations = 999)
centroids_df <- as.data.frame(fit_m$factors$centroids)
centroids_df$Level <- rownames(centroids_df)
ggplot(mushroomsc, aes(x = NMDS1, y = NMDS2)) +
geom_point(aes(color = class), alpha = 0.4, size = 1.5) +
geom_point(data = centroids_df, aes(x = NMDS1, y = NMDS2),
shape = 17, color = "black", size = 2.5,
inherit.aes = FALSE) +
geom_text_repel(data = centroids_df,
aes(x = NMDS1, y = NMDS2, label = Level),
size = 2.5, max.overlaps = 20,
inherit.aes = FALSE) +
scale_color_manual(values = c(e = "forestgreen", p = "firebrick"),
labels = c(e = "Edible", p = "Poisonous")) +
labs(title = "NMDS of Mushrooms with Factor Centroids",
color = "Class") +
coord_equal() +
theme_classic()
Not that readable….it might clear up w/o labels.
D)
We’ve got some filtering to do. Filter the centroids data frame to contain only variables with >0.5 \(R^2\) with respect to the ordination space. Then re-create your plot. If you want to eat a mushroom and live to tell about it, what characteristics should you look for? What characteristics should you avoid? (For full credit, you should spend some time on this response.)
r2 <- fit_m$factors$r
keep_vars <- names(r2)[r2 > 0.5]
sort(r2[keep_vars], decreasing = TRUE) # which variables made the cut ring.type stalk.surface.below.ring stalk.surface.above.ring
0.8167712 0.6858516 0.6804974
stalk.color.below.ring stalk.color.above.ring odor
0.5968820 0.5929885 0.5890816
spore.print.color
0.5084607 centroids_df$Variable <- fit_m$factors$var.id
centroids_filt <- centroids_df %>%
filter(Variable %in% keep_vars)
ggplot(mushroomsc, aes(x = NMDS1, y = NMDS2)) +
geom_point(aes(color = class), alpha = 0.4, size = 1.5) +
geom_point(data = centroids_filt, aes(x = NMDS1, y = NMDS2),
shape = 17, color = "black", size = 2.5,
inherit.aes = FALSE) +
geom_text_repel(data = centroids_filt,
aes(x = NMDS1, y = NMDS2, label = Level),
size = 3, max.overlaps = 30,
inherit.aes = FALSE) +
scale_color_manual(values = c(e = "forestgreen", p = "firebrick"),
labels = c(e = "Edible", p = "Poisonous")) +
labs(title = "NMDS of Mushrooms: Centroids for Variables with R² > 0.5",
color = "Class") +
coord_equal() +
theme_classic()
The plot shows that odor, spore print color, ring type, stalk surface, and stalk color help separate edible and poisonous mushrooms. The big poisonous cluster on the right is mostly associated with foul odor, chocolate spore prints, large rings, silky stalks, and buff or brown stalk colors. Another smaller poisonous group has a musty odor, no ring, and a cinnamon stalk. Edible mushrooms are more associated with almond or anise odors, smooth or fibrous stalks, gray stalk color, and certain ring types. However, some traits like no odor and white, black, or brown spore prints appear in both groups, so they aren’t as useful on their own.