Activity 2.1 - dissimilarity measures

setwd() SUBMISSION INSTRUCTIONS:

Submit both:

• Your source .qmd file;
• A link to the published version of your rendered html.

Question 1

A)

Consider the data set below:

Person	Age (years)	Income (k$)	Gender	Education
P1	25	35	Female	Bachelor
P2	40	50	Male	Bachelor
P3	30	65	Female	Master

Which distance measure would be best for measuring the dissimilarity between these individuals? Find the 3 pairwise distances using your selected metric “from scratch”. Which pair is most different? Most similar?

GOWER

\[ P1 - P2 = (15/15 + 15/30 + 1 + 0)/4 = .625\] \[ P1 - P3 = (5/15 + 30/30 + 0 + 1)/4 = .583\]

MOST SIMILAR

\[ P2 - P3 = (10/15 + 15/30 + 1 + 1)/4 = .792\]

MOST DIFFERENT

B)

Consider the data set below:

Person	Preferred Cuisine	Hobby	Device Used Most
P1	Italian	Reading	Laptop
P2	Mexican	Reading	Phone
P3	Italian	Hiking	Laptop

Which distance measure would be best for measuring the dissimilarity between these individuals? Find the 3 pairwise distances using your selected metric “from scratch”. Which pair is most different? Most similar?

JACCARD

	Italian	Mexican	Reading	Hiking	Laptop	Phone
P1	1	0	1	0	1	0
P2	0	1	1	0	0	1
P3	1	0	0	1	1	0

P1 P2	1	0
1	1	2
0	2	1

\[P1 - P2 = 4/5\]

P1 P3	1	0
1	2	1
0	1	2

\[P1 - P3 = 2/4\]

MOST SIMILAR

P2 P3	1	0
1	0	3
0	3	0

\[P2 - P3 = 6/6\] MOST DIFFERENT

C)

Consider the data set below representing counts of three different species at three locations:

Site	Species A	Species B	Species C
Site 1	4	2	0
Site 2	1	3	1
Site 3	0	5	2

Which distance measure would be best for measuring the dissimilarity between these individuals? Find the 3 pairwise distances using your selected metric “from scratch”. Which pair is most different? Most similar?

Bray-Curtis

\[1 and 2 = 3 + 1 + 1 = 5 / 11\]

\[1 and 3 = 4 + 3 + 2 = 9/13\]

MOST DIFFERENT

\[2 and 3 = 1 + 2 + 1 = 4/12\]

MOST DIFFERENT

Question 2

You’re a city manager in charge of building a new fire station. The fire station serves a city of only 4 households. There are 3 locations in the city suitable for the fire station.

The grid below shows the locations of the houses (grey squares) and suitable fire station locations (red circles):

A)

Which distance metric would be most appropriate for determining the optimal place to build the fire station? Explain your reasoning.

Manhattan. This is simulates the distance it will take to respond to emergencies via roads. The shorter the manhattan distance on average, the better, and in theory the shorter average response time.

B)

Which is the optimal place of the available locations to build the fire station that minimizes the cumulative distance between the station and the houses?

B

C)

Now consider a much larger city:

library(tidyverse)

large_city <- read.csv('Data/large_city.csv')


ggplot(large_city, aes(x = X, y = Y)) +
  geom_point(pch = 15) +           
  theme_bw(base_size = 14) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_line(color='lightgrey'),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  scale_x_continuous(expand = c(0.01,0.01), limits = c(0,100), breaks = 0:100) +
  scale_y_continuous(expand = c(0.01,0.01), limits = c(0,100), breaks = 0:100) 

The code below produces the coordinates of all locations in the city:

(possibilities <- expand.grid(X_FS = 0:100, Y_FS = 0:100) 
      %>% mutate(FireStationID = 1:n())
) %>% head

  X_FS Y_FS FireStationID
1    0    0             1
2    1    0             2
3    2    0             3
4    3    0             4
5    4    0             5
6    5    0             6

The next chunk of code crosses each house coordinate with all city coordinates:

distances <- large_city %>%
  crossing(possibilities) %>%
  mutate(dist = abs(X - X_FS) + abs(Y - Y_FS)) %>%
  group_by(FireStationID, X_FS, Y_FS) %>%
  summarize(total_dist = sum(dist), .groups = "drop") %>%
  anti_join(large_city, by = c("X_FS" = "X", "Y_FS" = "Y")) %>%
  filter(total_dist == min(total_dist))

distances_L2 <- large_city %>%
  crossing(possibilities) %>%
  mutate(dist = sqrt((X - X_FS)^2 + (Y - Y_FS)^2)) %>%
  group_by(FireStationID, X_FS, Y_FS) %>%
  summarize(total_dist = sum(dist), .groups = "drop") %>%
  anti_join(large_city, by = c("X_FS" = "X", "Y_FS" = "Y")) %>%
  filter(total_dist == min(total_dist))

ggplot(large_city, aes(x = X, y = Y)) +
  geom_point(pch = 15) +
  geom_point(data = distances, aes(x = X_FS, y = Y_FS), color = "red", size = 1) +
  geom_point(data = distances_L2, aes(x = X_FS, y = Y_FS), color = "blue", size = 1) +
  theme_bw(base_size = 14)

Suppose any location where a house is not already built is a candidate for a fire station location. Complete the dplyr chain above to find the fire station location that minimizes the L1 distance from all houses (there may be multiple optimal locations!). Add the location(s) of the fire station to the graph of the city using red dot(s). You should of course make sure that you don’t build a fire station on top of an existing house!

D)

Modify your code to find the optimal location of the station if L2 distance is used instead. Add the new location to the grid.

Question 3

Consider the R code below which produced one of the plots from the lectures slides:

library(plotly)

# Read and wrangle
city77 <- (read.csv("Data/City77.csv")
        %>% mutate(City = gsub('\\.',' ',City))
        %>% mutate(highlight = ifelse(City %in% c("Minneapolis MN", "Kansas City MO", "Milwaukee WI"), 'yes','no'))
) 

#Plot
plot_ly(city77,
        x = ~popdens,
        y = ~k12enr,
        z = ~medinc,
        type = "scatter3d",
        mode = "markers+text",
        text = ~ifelse(highlight == "yes", City, ""),
        textposition = "top center",
        color = ~highlight,
        colors = c("gray", "red"),
        marker = list(size = 5)) %>%
  layout(title = "3D Plot: Population Density vs K-12 vs Median Income",
         scene = list(
           xaxis = list(title = "Population Density"),
           yaxis = list(title = "K-12 enrollment %"),
           zaxis = list(title = "Median Income")
         ))

Use this code as a template to analyze bankruptcy data (source: Johnson and Wichern 6e page 657). This data set contains four metrics on both bankrupt and financially sound firms. The metrics were recorded 2 years prior to bankruptcy for the bankrupt firms, and about the same time for the financially sound firms.

Variables:

• bankrupt = 1 for yes, 0 for financially sound
• x1 = cash flow / total debt
• x2 = net income/total assets
• x3 = current assets/current liabilities
• x4 = current assets/net sales

bankruptcy <- read.csv('Data/bankruptcy_data.csv')

A)

Create an interactive 3D scatterplot of this data using x1, x2, and x3. Color-code the points by bankruptcy status. Imagine your boss is a bank manager. Use your plot to explain what are some warning signs to look out for in terms of the characteristics of bankrupt banks.

plot_ly(bankruptcy,
        x = ~x1,
        y = ~x2,
        z = ~x3,
        type = "scatter3d",
        mode = "markers",
        color = ~factor(bankrupt),
        colors = c("gray", "red"),
        marker = list(size = 5)) %>%
  plotly::layout(title = "3D Plot: Cash Flow/TD vs Net Income/TA vs Current Ratio",
                 scene = list(
                   xaxis = list(title = "Cash Flow / Total Debt"),
                   yaxis = list(title = "Net Income / Total Assets"),
                   zaxis = list(title = "Current Assets / Current Liabilities")
                 ))

The chart shows a clear warning sign when there is low (or negative) cash flow and low current assets.

B)

Find the mean vectors of x1 through x4 for bankrupt and sound banks (so 2 mean vectors). Then consider a new bank with the following metrics: x1 = 0.31, x2 = 0.04, x3 = 4.47, x4 = 0.30

Is this bank doing well or are they in trouble? Justify your answer using distance between this bank’s metrics and the mean vector for each bank type.

mean_vectors <- bankruptcy %>%
  group_by(bankrupt) %>%
  summarize(across(c(x1, x2, x3, x4), mean), .groups = "drop")

mean_vectors

# A tibble: 2 × 5
  bankrupt      x1      x2    x3    x4
     <int>   <dbl>   <dbl> <dbl> <dbl>
1        0  0.235   0.0556  2.59 0.427
2        1 -0.0690 -0.0814  1.37 0.438

new_bank <- c(x1 = 0.31, x2 = 0.04, x3 = 4.47, x4 = 0.30)

sound_mean <- mean_vectors %>% filter(bankrupt == 0) %>% select(x1:x4) %>% as.numeric()
bankrupt_mean <- mean_vectors %>% filter(bankrupt == 1) %>% select(x1:x4) %>% as.numeric()

dist_sound <- sqrt(sum((new_bank - sound_mean)^2))
dist_bankrupt <- sqrt(sum((new_bank - bankrupt_mean)^2))

dist_sound

[1] 1.882231

dist_bankrupt

[1] 3.131779

The new banks metrics are closer to the mean vector of banks of financially sound banks (1.88) than bankrupt banks (3.13) which suggests the bank is currently in good financial standing.

Question 4

These data come from a collaboration with WSU biology faculty and students and published in Water. The Whitewater River has three forks: north, middle, and south.

Image source: Water

Electrofishing at over 60 locations on the three forks was used to sample fish populations. The data are below.

whitewater <- read.csv('Data/Whitewater.csv')

A)

Explain why Bray-Curtis is appropriate for measuring similarity/dissimilarity in populations across sites. Use Bray-Curtis and data management techniques to find the two most similar, and the two most dissimilar sites.

Bray-Curtis ignores double-zeros which helps us out a lot here becasue many of the species are absent which will create double-zeros.

library(vegan)

Loading required package: permute

library(tidyverse)

species_data <- whitewater %>%
  select(-fork, -site, -siteID, -lat, -long)
bc_dist <- vegdist(species_data, method = "bray")

bc_matrix <- as.matrix(bc_dist)

bc_long <- as.data.frame(as.table(bc_matrix)) %>%
  rename(site1 = Var1, site2 = Var2, dissimilarity = Freq) %>%
  filter(site1 != site2) %>%
  filter(as.integer(site1) < as.integer(site2))

bc_long %>% slice_min(dissimilarity, n = 2)

  site1 site2 dissimilarity
1    23    24    0.08038585
2    46    47    0.12720848

bc_long %>% slice_max(dissimilarity, n = 2)

   site1 site2 dissimilarity
1      1     2             1
2      1     6             1
3      1     8             1
4      1    10             1
5      1    11             1
6      1    13             1
7      1    15             1
8      1    17             1
9      1    18             1
10     3    18             1
11     1    19             1
12     1    20             1
13     3    20             1
14     1    21             1
15     3    21             1
16     1    22             1
17     2    22             1
18     3    22             1
19     4    22             1
20     7    22             1
21     1    23             1
22     2    23             1
23     3    23             1
24     4    23             1
25     7    23             1
26     1    24             1
27     2    24             1
28     3    24             1
29     4    24             1
30     7    24             1
31    18    25             1
32    20    25             1
33    21    25             1
34    22    25             1
35    23    25             1
36    24    25             1
37     1    32             1
38     1    33             1
39     1    34             1
40     3    34             1
41    25    34             1
42     1    35             1
43     2    35             1
44     3    35             1
45    25    35             1
46     1    36             1
47     3    36             1
48    25    36             1
49     1    37             1
50     1    38             1
51     3    38             1
52    25    38             1
53     1    39             1
54     2    39             1
55     3    39             1
56     4    39             1
57    25    39             1
58     1    40             1
59     3    40             1
60    25    40             1
61    21    41             1
62    22    41             1
63    23    41             1
64    24    41             1
65    32    41             1
66    33    41             1
67    34    41             1
68    35    41             1
69    36    41             1
70    37    41             1
71    38    41             1
72    39    41             1
73    22    42             1
74    23    42             1
75    24    42             1
76    22    43             1
77    23    43             1
78    24    43             1
79    22    44             1
80    23    44             1
81    24    44             1
82     1    47             1
83     1    48             1
84     1    49             1
85     1    51             1
86     1    52             1
87     1    53             1
88     1    54             1
89     1    55             1
90     1    56             1
91     1    57             1
92     1    59             1
93     3    59             1
94    25    59             1
95     1    60             1
96     1    61             1

Most similar:

Most dissimlar: sites 23 & 24 and 46 & 47

Most similar: 96 different pairs exactly the same

B)

The code below transforms the data set to measure only presence/absence of each species:

present <- \(species) ifelse(species>0, 1, 0)

whitewater_present <- (whitewater
                       %>% mutate(across(.cols=American.brook.lamprey:slimy.sculpin, 
                          .fns=present))
                      ) 

Consider these two sites:

whitewater_present %>% 
  slice(2, 4)

  fork            site siteID      lat      long American.brook.lamprey
1   NF  Pries source 2   NF33 44.04202 -92.26845                      0
2   NF         Pries 2   NF31 44.04505 -92.26420                      0
  Rainbow.trout Brown.trout Brook.trout Central.stoneroller Common.shiner
1             0           0           0                   0             0
2             0           0           0                   0             0
  Sand.shiner Southern.redbelly.dace Bluntnose.minnow Fathead.minnow
1           0                      0                0              0
2           0                      1                0              0
  Blacknose.dace Longnose.dace Creek.chub White.sucker Brook.stickleback
1              1             0          1            0                 0
2              1             1          1            0                 1
  Green.sunfish Bluegill Fantail.darter Johnny.darter mottled.sculpin
1             0        0              1             1               0
2             0        0              1             1               0
  slimy.sculpin
1             0
2             0

Compute the Jaccard dissimilarity metric between these two sites from scratch.

	1	0
1	4	0
0	3	14

3/7 = .429

C)

Use the Gower measure on the count data, and also find the Jaccard dissimilarities for all sites using the present/absent data. Turn each into a data frame and join them with the Bray-Curtis measures (so you have 3 measures for each unique pair of site comparisons). Then create pairwise scatterplots of the 3 dissimilarity measures versus each other. How do they compare?

library(vegan)
library(tidyverse)
library(patchwork)


bc_dist <- vegdist(species_data, method = "bray")
bc_matrix <- as.matrix(bc_dist)
rownames(bc_matrix) <- colnames(bc_matrix) <- whitewater$site

gower_matrix <- as.matrix(vegdist(species_data, method = "gower"))
rownames(gower_matrix) <- colnames(gower_matrix) <- whitewater$site

present_data <- whitewater_present %>% select(American.brook.lamprey:slimy.sculpin)
jaccard_matrix <- as.matrix(vegdist(present_data, method = "jaccard", binary = TRUE))
rownames(jaccard_matrix) <- colnames(jaccard_matrix) <- whitewater$site

matrix_to_long <- function(mat, colname) {
  as.data.frame(as.table(mat)) %>%
    rename(site1 = Var1, site2 = Var2, !!colname := Freq) %>%
    filter(as.character(site1) < as.character(site2))
}

all_measures <- matrix_to_long(bc_matrix, "bray_curtis") %>%
  inner_join(matrix_to_long(gower_matrix, "gower"), by = c("site1", "site2")) %>%
  inner_join(matrix_to_long(jaccard_matrix, "jaccard"), by = c("site1", "site2"))

p1 <- ggplot(all_measures, aes(bray_curtis, jaccard)) +
  geom_point(alpha = 0.5) +
  labs(title = "Bray-Curtis vs Jaccard") +
  theme_bw()

p2 <- ggplot(all_measures, aes(bray_curtis, gower)) +
  geom_point(alpha = 0.5) +
  labs(title = "Bray-Curtis vs Gower") +
  theme_bw()

p3 <- ggplot(all_measures, aes(gower, jaccard)) +
  geom_point(alpha = 0.5) +
  labs(title = "Gower vs Jaccard") +
  theme_bw()

p1 / p2 / p3

Bray-Curtis vs Jaccard shows the most agreement, Bray-Curtis vs Gower shows some agreement but there are problems when bray-curtis = 1.00, and Gower vs Jaccard is not close.

D)

Find the 2 pairs of sites with the largest discrepancies between the Bray-Curtis and Jaccard measures. Can you identify why they are in disagreement?

discrepancies <- all_measures %>%
  mutate(diff = abs(bray_curtis - jaccard)) %>%
  arrange(desc(diff))

discrepancies %>%
  select(site1, site2, bray_curtis, jaccard, diff) %>%
  slice(1:2)

         site1          site2 bray_curtis jaccard      diff
1       Borgen        BReiter   0.8014184       0 0.8014184
2  Irke Source  Shea Crossing   0.7225951       0 0.7225951

Both pairs have a Jaccard dissimilarity of 0, meaning the sites have the exact same species. However, both pairs have high Bray-Curtis dissimilarity, meaning the species are present in different amounts at each site. This shows that Jaccard only considers whether species are present, while Bray-Curtis also considers how abundant each species is.