college <- read.csv('Data/College.csv')Activity 2.2 - multidimensional scaling
SUBMISSION INSTRUCTIONS
Submit both:
- Your .qmd;
- A link to your rendered html.
Question 1
Suppose we have 3 data points in 2 dimensions:
A)
Using taxicab as the distance metric, find all pairwise distances between points in the full dimension space.
- \(d(A,B) = 4\); \(d(A,C) = 6\), \(d(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?
- First ordination:
- Distances: \(d(A,B) = 3\); \(d(A, C) = 4\); \(d(B, C) = 1\)
- Stress: \(\sqrt{(4-3)^2 + (6-4)^2 + (4-1)^2} = 3.74\)
- Second ordination:
- Distances: \(d(A,B) = 1\); \(d(A, C) = 2\); \(d(B, C) = 3\)
- Stress: \(\sqrt{(4-1)^2 + (6-2)^2 + (4-3)^2} = 5.1\)
- Third ordination:
- Distances: \(d(A,B) = 2\); \(d(A, C) = 4\); \(d(B, C) = 6\)
- Stress: \(\sqrt{(4-2)^2 + (6-4)^2 + (4-6)^2} = 3.46\)
So the third ordination is the best 1D representation of these 2D data, as it has the lowest stress.
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?
Properties of this plot:
- Variables A and C are quite strongly (and positively) correlated with each other, as they point in similar directions.
- Variables A and C are strongly negatively correlated with B, as these arrows point in opposite directions..
- A, C and B all have weak correlations with D; in fact it looks like \(Cor(B,D)\) and \(Cor(A,D)\) might be near 0 due to the orthogonality (right angles) of these arrows.
## 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.00We can rule out Matrix 1 right away because it has a positive correlation between B and C, while we know B and C are negatively correlated.
## 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.00Matrix 2 can be ruled out since \(Cor(B,D) = 0.39\) while \(Cor(C,D) = -0.06\): the biplot certainly is not consistent with B and D be more correlated than C and D.
## 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.00We can rule out Matrix 3 right away due to \(Cor(A,C) = -0.94\).
## 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 is the one! We have:
- perfect negative correlation between A and B (those arrows point in opposite directions);
- strong positive correlation between A and C (those arrows point in very similar directions);
- Near 0 correlation between A and D (those arrows are almost at a right angle);
- Near 0 correlation between B and D (those arrows are almost at a right angle);
- Strong negative correlation between B and C (these arrows point in almost opposite directions);
- Weak positive correlation between C and D.
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.00- A and B are weakly correlated; these arrows should be pointing at near right angles.
- A and C are weakly correlated these arrows should be pointing at near right angles.
- B and C are strongly correlated these arrows should be pointing in similar directions.
Plot 4 is most consistent with this correlation matrix.
Now consider four biplots below showing the 3-dimensional data set reduced to 2 dimensions using NMDS:
Which plot does the correlation matrix belong to?
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.
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?
Euclidean or Manhattan distance is appropriate as we have numeric variables that aren’t counts. We’ll use Euclidean, which does require scaling to prevent large-scale variables from dominating the distance calculations.
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.
library(vegan)
library(tidyverse)
(numeric_scaled_vars <- college
%>% select(Apps:Grad.Rate)
%>% mutate(across(everything(), scale))
) %>% head Apps Accept Enroll Top10perc Top25perc F.Undergrad
1 -0.3466585 -0.32099869 -0.06346802 -0.2584163 -0.1917039 -0.1680076
2 -0.2107483 -0.03867808 -0.28839845 -0.6552335 -1.3530399 -0.2096534
3 -0.4066037 -0.37607569 -0.47781355 -0.3151045 -0.2926897 -0.5492116
4 -0.6678304 -0.68124306 -0.69198176 1.8390460 1.6765321 -0.6576558
5 -0.7257086 -0.76406254 -0.78023197 -0.6552335 -0.5956469 -0.7114656
6 -0.6239051 -0.62820595 -0.66938109 0.5919062 0.3132247 -0.6230197
P.Undergrad Outstate Room.Board Books Personal PhD Terminal
1 -0.20907246 -0.7458755 -0.9642836 -0.6019244 1.2692276 -0.1629230 -0.1156542
2 0.24414978 0.4572019 1.9079789 1.2150971 0.2353633 -2.6739232 -3.3760014
3 -0.49677006 0.2011751 -0.5539604 -0.9047614 -0.2594146 -1.2040694 -0.9307410
4 -0.52041644 0.6262293 0.9961495 -0.6019244 -0.6877298 1.1844430 1.1748999
5 0.00899969 -0.7160471 -0.2165835 1.5179340 0.2353633 0.2045405 -0.5231976
6 -0.53486701 0.7604569 -0.9323696 -0.2990875 -0.9831196 -0.3466547 -0.4552737
S.F.Ratio perc.alumni Expend Grad.Rate
1 1.0131234 -0.8670157 -0.50158700 -0.3180471
2 -0.4773970 -0.5442217 0.16600292 -0.5509070
3 -0.3005556 0.5855576 -0.17717583 -0.6673369
4 -1.6142346 1.1504472 1.79169738 -0.3762621
5 -0.5531862 -1.6740009 0.24164779 -2.9377211
6 -1.1847626 -0.9477142 0.01279812 -0.6091220euc_dist <- vegdist(numeric_scaled_vars, method='euclidean')
college_mds <- metaMDS(euc_dist, k = 2, trymax = 20, trace = 0)
college_update <- bind_cols(college, college_mds$points)
college_mds$stress[1] 0.1285852The stress of this solution is 0.13, which is acceptable. That means we are getting a decent representation of the similarities and dissimilarities between these colleges with respect to all the numeric variables with a much simpler 2-dimensional representation.
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?
low_dimensional_distances <- vegdist(college_mds$points, method='euclidean')
stressdf <- data.frame(fullD_distance = as.vector(euc_dist),
lowD_distance = as.vector(low_dimensional_distances)
)
cor(stressdf) fullD_distance lowD_distance
fullD_distance 1.0000000 0.9638556
lowD_distance 0.9638556 1.0000000ggplot(data = stressdf, aes(x = fullD_distance, y = lowD_distance))+
geom_point(shape='.')+
labs(x = 'Full dimensional dissimilarity',
y= 'Low dimensional dissimilarity',
title = 'Stress plot of the 2D configuration',
subtitle = 'Correlation = 0.96'
)+
theme_classic()
cor(stressdf) fullD_distance lowD_distance
fullD_distance 1.0000000 0.9638556
lowD_distance 0.9638556 1.0000000The correlation between the full-dimensional and 2-dimensional dissimilarities is 0.96, which is quite high. This is a good thing; we want the 2-dimensional dissimilarities to be a faithful representation of the full-dimensional distances, that’s the point of NMDS!
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?
college_fit <- envfit(ord = college_mds, env = numeric_scaled_vars)
arrowdf <- college_fit$vectors$arrows %>%
data.frame %>%
rownames_to_column('variable') %>%
mutate(r2 = college_fit$vectors$r) %>%
mutate(NMDS1 = NMDS1*sqrt(r2)*4,NMDS2 = NMDS2*sqrt(r2)*4)
arrowdf %>% slice_max(n=3, order_by = r2) variable NMDS1 NMDS2 r2
1 F.Undergrad 3.133711 -1.3932492 0.7350805
2 Accept 3.348509 -0.7017151 0.7315573
3 Apps 3.417185 -0.1311865 0.7308978arrowdf %>% slice_min(n=3, order_by = r2) variable NMDS1 NMDS2 r2
1 Books 0.1290485 0.2599854 0.00526537
2 Personal 0.8971931 -1.5298296 0.19658339
3 Room.Board 0.6416303 2.3488534 0.37055010Most responsible for separating universities:
- number of full time undergraduates;
- instructional expenditure per student;
- Number of applications received and accepted, and number of new students enrolled.
Least responsible for separating universities:
- estimated book costs,
- estimated personal spending,
- Room and board costs
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?
Plot below, note that I am filtering to variables that have at least a 60% \(R^2\) with the ordination space:
important_arrows <- arrowdf %>% filter(r2 > 0.6)
library(ggrepel)
ggplot(data = college_update, aes(x = MDS1, y = MDS2)) +
geom_point(aes( color = Private), alpha = 0.5) +
stat_ellipse(geom='polygon', aes(fill = Private), alpha = 0.4, level = 0.95)+
geom_segment(aes(x = 0, y = 0,
xend = NMDS1, yend = NMDS2),
arrow = arrow(length = unit(.2, 'cm')),
data = important_arrows) +
geom_text_repel(aes(x = NMDS1, y = NMDS2, label = variable),
data = important_arrows) +
theme_classic(base_size = 14) 
It appears that number of undergraduates, applications, accepted, and enrolled students (basically, the size of the university) separates points along the horizontal dimension. The vertical dimension seems to be relying on the cost and selectivity of the university: out-of-state tuition, expenditure per student, and percent of new students from the top 10% or 25% of class is what’s represented by the vertical dimension.
To summarize this plot in a nutshell, the public universities tend to be large, while the private universities tend to be expensive and selective.
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.
The first step is to pivot, then use the pivoted data to created a faceted boxplot:
(college_long <- college
%>% pivot_longer(cols = Apps:Grad.Rate,
values_to = 'value',
names_to='variable'
)
%>% inner_join(important_arrows, by = 'variable')
%>% mutate(variable = fct_reorder(variable, r2, .desc=T))
)%>% head# A tibble: 6 × 7
Name Private variable value NMDS1 NMDS2 r2
<chr> <chr> <fct> <dbl> <dbl> <dbl> <dbl>
1 Abilene Christian University Yes Apps 1660 3.42 -0.131 0.731
2 Abilene Christian University Yes Accept 1232 3.35 -0.702 0.732
3 Abilene Christian University Yes Enroll 721 3.22 -1.14 0.728
4 Abilene Christian University Yes Top10perc 23 1.63 2.86 0.678
5 Abilene Christian University Yes Top25perc 52 1.79 2.53 0.602
6 Abilene Christian University Yes F.Undergrad 2885 3.13 -1.39 0.735ggplot(data = college_long) +
geom_boxplot(aes(x = Private, y = value, fill=Private)) +
facet_wrap(~variable, scales = 'free_y')+
theme_classic()
As we can see, public universities tend to have higher number of undergraduates, applications, accepted, and enrolled students; private universities have higher expenditure, higher percentage of students at the top of their class, and higher out-of-state tuition and expenditure per student.
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.
library(fastDummies)
mushrooms_onehot <- mushroomsc %>%
select(-class) %>%
dummy_cols(., remove_selected_columns = TRUE)
jaccard_dist <- vegdist(mushrooms_onehot, method='jaccard')B)
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?
mush_mds2 <- metaMDS(jaccard_dist, k = 2, trymax=1)Run 0 stress 0.1459574
Run 1 stress 0.148828
*** Best solution was not repeated -- monoMDS stopping criteria:
1: scale factor of the gradient < sfgrminmush_mds3 <- metaMDS(jaccard_dist, k = 3, trymax=1)Run 0 stress 0.09984206
Run 1 stress 0.09992898
... Procrustes: rmse 0.005100831 max resid 0.05699259
*** Best solution was not repeated -- monoMDS stopping criteria:
1: scale factor of the gradient < sfgrminmush_mds2$stress[1] 0.1459574mush_mds3$stress[1] 0.09984206The 3D configuration has a stress of 0.0998, in the “very good” realm, while the 2D is 0.146, merely “acceptable.”
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?
#update the mushrooms data frame with MDS coordinates:
mushrooms_with_mds <- bind_cols(mushroomsc, mush_mds2$points)
#Fit mushroom characteristics on ordination:
mushrooms_noclass <- mushroomsc %>% select(-class)
mush_fit <- envfit(ord= mush_mds2, env = mushrooms_noclass)
#Extract centroids to data frame:
centroids <- mush_fit$factors$centroids %>%
data.frame %>%
rownames_to_column('variable')
#Produce plot:
library(ggrepel)
ggplot(mushrooms_with_mds) +
geom_point(aes(x = MDS1, y = MDS2, color = class),size = 3, alpha = 0.7) +
geom_point(aes(x = NMDS1, y = NMDS2), data = centroids, shape =15 ) +
geom_text_repel(aes(x = NMDS1, y = NMDS2, label = variable), data = centroids)+
labs(x = "MDS Dimension 1", y = "MDS Dimension 2") +
theme_classic(base_size = 14)
There are way too many centroid points on this plot to make sense of it!
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.)
#Sort the variables by R2:
mush_fit$factors$r %>%
sort(., decreasing = TRUE) ring.type stalk.surface.above.ring stalk.surface.below.ring
0.81858434 0.68658045 0.67062433
stalk.color.below.ring stalk.color.above.ring odor
0.60082927 0.59836672 0.58485776
spore.print.color stalk.root population
0.53025223 0.40341466 0.39448550
bruises habitat gill.color
0.34218841 0.31598609 0.30981461
cap.color gill.spacing cap.surface
0.26634394 0.26255209 0.16811246
cap.shape gill.size stalk.shape
0.07578453 0.07551151 0.06948195
ring.number
0.03935608 #Filter to the ones with R2 > 0.5, which is the first 7:
important <- mush_fit$factors$r %>%
sort(., decreasing = TRUE)%>%
head(7)
# Now filter the centroids data.
centroids_important <- centroids %>%
#Mutate to get rid of the level indicator at end of string:
mutate(variable_only = str_sub(variable, start=1, end = -2)) %>%
# Filter to variables that are in important vector:
filter(variable_only %in% names(important))ggplot(mushrooms_with_mds) +
geom_point(aes(x = MDS1, y = MDS2, color = class),size = 3, alpha = 0.7) +
geom_point(aes(x = NMDS1, y = NMDS2), data = centroids_important, shape =15 ) +
geom_text_repel(aes(x = NMDS1, y = NMDS2, label = variable), data = centroids_important,
max.overlaps = 20)+
labs(x = "MDS Dimension 1", y = "MDS Dimension 2") +
theme_classic(base_size = 14)
If you want to eat a mushroom and not die, you should definitely avoid yellow, buff, or brown stalks both above and below the ring. Gray and pink is much safer! Also avoid foul, pungent, or musty-smelling mushrooms, or odorless mushrooms and go for the ones that smell like anise or almond! Don’t be fooled by the silky stalk surfaces: they’re deadly. Fibrous, smooth stalk surfaces are your friend.