---
title: "💻 LSE 202A 2023: Week 04 - Lab"
author: "MY NAME (MY CANDIDATE NUMBER)"
date: 27 November 2023
format: html
self-contained: true
---

GENIAL Participants:

-   **Fill out this brief survey at the end of the lab:** 🔗 [link](https://forms.office.com/e/h0dXriciyy) (requires LSE login)
  - You don't have to use ChatGPT today, but if you do, please share the chat log!

# ⚙️ Setup

Load the libraries for today's lab.

```{r, message=FALSE, warning=FALSE, echo=FALSE}
# Tidyverse packages we will use
library(ggplot2)
library(dplyr)     
library(tidyr)     
library(readr)     
library(forcats)
library(stringr)

library(recipes)
library(broom)

# Non-tidymodels packages
library(dbscan)
```

## Part I - Meet a new dataset (15+5 mins)

No need to wait! Start tackling the action points below when you come to the classroom.

## 🎯 ACTION POINTS (15 min)

1. Download the data set:

```{r, message=FALSE, warning=FALSE}
# Download the zipped file
download.file("https://github.com/rudeboybert/JSE_OkCupid/raw/master/profiles_revised.csv.zip",
              dest="profiles_revised.csv.zip")

# Unzip it so we can use a profiles_revised.csv file
unzip("profiles_revised.csv.zip")
```

Keep a tab open with [data codebook](https://github.com/rudeboybert/JSE_OkCupid/blob/master/okcupid_codebook_revised.txt).

2. Read it and take a look:

```{r, message=FALSE, warning=FALSE}
df_okcupid <- 
  read_csv("profiles_revised.csv") %>%
  # Fix issues with the data, converting "does&rsquo;t" to "doesnt"
  mutate(across(where(is.character), ~ str_replace_all(., stringr::fixed("&rsquo;"), ""))) %>%
  mutate(across(where(is.character), ~ str_replace_all(., stringr::fixed("++"), "pp"))) %>%
  mutate(zodiac_sign=purrr::map_chr(sign, ~ str_split(., " ")[[1]][1]),
         zodiac_importance=case_when(str_detect(sign, "but it doesnt matter") ~ "low importance",
                                     str_detect(sign, "and its fun to think about") ~ "medium importance",
                                     str_detect(sign, "matters a lot") ~ "high importance"),
         religion_importance=case_when(str_detect(religion, "and laughing about it") ~ "very low importance",
                                       str_detect(religion, "but not too serious about it") ~ "low importance",
                                       str_detect(religion, "and somewhat serious about it") ~ "medium importance",
                                       str_detect(religion, "and very serious about it") ~ "high importance"),
         religion=purrr::map_chr(religion, ~ str_split(., " ")[[1]][1])) %>%
  mutate(income=if_else(income == -1, "not informed", sprintf("%d", income))) %>%
  select(-c(sign))

df_okcupid %>% glimpse()
```

3. **Refine and Tidy Up**: A closer look at the dataset reveals that certain columns, like `ethnicity`, `offspring`, and `speaks`, contain values that reflect multiple selected responses. This occurs because respondents were given the option to choose more than one answer:

```{r, message=FALSE, warning=FALSE}
df_okcupid %>% select(ethnicity, offspring, speaks) %>% head() %>% knitr::kable()
```

We want to convert those to **dummy variables** to better use them later. We cannot use recipes for several reasons (ask Jon if you want to understand them). Instead, we must employ advanced data-reshaping techniques to create dummies of these particular columns.

⚠️ **The code below is very advanced. The kind of which only the most advanced students of our sister course about data manipulation, DS105, would be able to understand (provided they know R)! So don't worry about understanding it all that well.**


```{r}
df_okcupid <-
  df_okcupid %>% 
  # After splitting the string, pivot to longer format (aka. melt the data frame) 
  separate_longer_delim(ethnicity, ", ") %>% 
  separate_longer_delim(offspring, ", ") %>%
  separate_longer_delim(speaks, ", ") %>%
  mutate(temporary_var=TRUE) %>%
  pivot_longer(c(ethnicity, offspring, speaks), names_to="multiple_choice_var", values_to="response") %>%
  distinct() %>%
  drop_na(response) %>%
  mutate(response=str_replace_all(response, stringr::fixed("'"), ""),
         response=str_replace(response, stringr::fixed(" ("), "_"),
         response=str_replace(response, stringr::fixed(")"), "")) %>%
  mutate(across(where(is.character), ~ str_replace_all(., stringr::fixed(" / "), "_or_"))) %>%
  mutate(across(where(is.character), ~ str_replace_all(., " ", "_"))) %>%
  pivot_wider(names_from = c(multiple_choice_var, response), 
              values_from = temporary_var,
              names_sep = "_",
              values_fill = FALSE)
```

What we want you to appreciate is that instead of a single column, `ethnicity`, we have multiple ethnicity-related ones:

```{r}
df_okcupid %>% select(starts_with("ethnicity")) %>% head()
```

The same applies to `speaks_` and `offspring_` variables.

4. Investigate the character columns. Just how many unique values are there in each of them?

```{r}
apply(df_okcupid %>% select(where(is.character)), 2, function(x){length(unique(x))})
```

Let's convert them all to factors!

```{r}
df_okcupid <- df_okcupid %>% mutate(across(where(is.character), as.factor))
```

5. **More dummy variables!** Let's use a recipe to convert the remaining categorical columns to dummy variables:

```{r}
dummy_recipe <-
  # We don't have an outcome variable here
  recipe(~., data = df_okcupid) %>%
  # Remove variables with more than 50% missing values
  step_filter_missing(all_predictors(), threshold=0.5) %>%
  # Remove rows with missing values
  step_naomit(all_predictors(), skip=FALSE) %>%
  # Remove zero variance predictors as they don't carry any information
  step_zv(all_numeric_predictors()) %>%
  # Normalize the data to avoid out-sized impact of some variables
  step_normalize(all_numeric_predictors()) %>%
  step_dummy(all_factor_predictors()) %>%
  # The dummy variables we created for ethnicity, offspring and speaks must be transformed to integers
  step_mutate(across(where(is.logical), as.integer)) %>%
  prep()

dummy_recipe %>% tidy()
```

5. Feel overwhelmed by taking a look at the resulting 458 columns!

```{r}
summary(dummy_recipe) %>% arrange(variable)
```

### 🗣CLASSROOM DISCUSSION (5 min)

(Your class teacher will mediate this discussion)

- How would you explore this data if dimensionality reduction was not an option?
- How would you answer: what are the **most common types of profiles** one can find on OK Cupid?

> Write your answers and reflections here

## Part II - The power of PCA (30 min)

### 🎯 ACTION POINTS

1. To try to make sense of the sheer number of combinations of attributes, we will run PCA:

```{r}
pca_recipe <-
  # Reuse existing recipe
  dummy_recipe %>%
  # I'm not sure yet how many components I want to keep; let's say 30.
  step_pca(all_predictors(), num_comp = 30, keep_original_cols=TRUE) %>%
  prep()


# Apply the recipe to our data set
pca_results <- bake(pca_recipe, df_okcupid)
```

2. How much 'information' are we keeping after compressing the data with PCA? 

```{r out.width=7, out.width=4, echo=FALSE}
pca_step <- pca_recipe %>% tidy() %>% filter(type == "pca") %>% pull(number)

pca_fit <- pca_recipe$steps[[pca_step]]$res

# the eigenvalues determine the amount of variance explained
eigenvalues = 
  pca_fit %>%
  tidy(matrix = "eigenvalues")

eigenvalues %>%
  filter(PC <= 5) %>%
  ggplot(aes(PC, cumulative)) +
  geom_col(fill="#C63C4A") +
  geom_text(aes(label=scales::percent(cumulative, accuracy=0.01)),
            nudge_y=0.05, fontface="bold", size=rel(4.5)) + 
  
  # Customise scales
  scale_y_continuous(name="Cumulative Explained Variance", 
                     labels=scales::percent,
                     limits=c(0, 1), 
                     breaks=seq(0, 1, 0.2)) +
  scale_x_continuous(name="Principal Component") +
  
  labs(title="The first PC explains ~33% of the variance in the data. Good compression!",
       subtitle="The second PC adds some 7% more, and so on...",
       caption="Figure 1. Explained variance\nper component") +
  
  # Prettify plot a bit
  theme_bw() +
  theme(plot.title=element_text(size=rel(1.2)),
        plot.subtitle = element_text(size=rel(1)),
        axis.title=element_text(size=rel(1.3)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        axis.text=element_text(size=rel(1.2)))
```

3. Let's focus on the first two components, as they are common plotting practices.

Think of the plot below as a ~37% compressed 2D representation of our 400+ numeric/dummy variables.

```{r out.width=7, out.width=4, echo=FALSE}
ggplot(pca_results, aes(PC01, PC02)) +
  geom_point(alpha=0.2, size=rel(2.5)) +
  labs(title = "There is a huge blob in the middle! I wonder what that represents.", 
       subtitle = "There are some clear outliers here",
       caption = "Figure 2. Distribution of samples using\nPC01 and PC02 as coordinates.") +
  
  scale_y_continuous(breaks=seq(-7.5, 7.5, 2.5), limits=c(-7.5, 7.5)) +
  
    # Prettify plot a bit
  theme_bw() +
  theme(plot.title=element_text(size=rel(1.2)),
        plot.subtitle = element_text(size=rel(1)),
        axis.title=element_text(size=rel(1.3)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        axis.text=element_text(size=rel(1.2)))
```

4. How concentrated is this main blob?

Let's change the `geom()` to see density instead of the dots:

```{r out.width=10, out.width=4, echo=FALSE}
ggplot(pca_results, aes(PC01, PC02)) +
  geom_bin2d() +
  labs(title = "Here is how concentrated the main blob is", 
       subtitle = "It's time to put your understanding of plots to the test!",
       caption = "Figure 3. Density of samples using\nPC01 and PC02 coordinates") +
  
  scale_fill_viridis_b() +
  # I tweaked those manually after seeing the plot
  scale_y_continuous(breaks=seq(-7.5, 7.5, 2.5), limits=c(-7.5, 7.5)) +
  
  # Prettify plot a bit
  theme_bw() +
  theme(plot.title=element_text(size=rel(1.2)),
        plot.subtitle = element_text(size=rel(1)),
        axis.title=element_text(size=rel(1.3)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        axis.text=element_text(size=rel(1.2))) 
```

### 🗣CLASSROOM DISCUSSION

(Your class teacher will mediate this discussion)

- How could you use the plots above to investigate the most common types of OK Cupid profiles? 

> Write your answers and reflections here

### 🧑‍🏫 TEACHING MOMENT

Your class teacher will now guide the conversation and explain the plot below. If needed, they will recap how PCA works.

```{r fig.width=10, out.width=8, echo=FALSE}
tidied_pca <- tidy(pca_recipe, pca_step)

plot_df <- 
  tidied_pca %>%
  filter(component %in% paste0("PC", 1:5)) %>%
  arrange(desc(abs(value))) %>% 
  group_by(component) %>% 
  # Get 10 biggest contributors to this PC
  slice_head(n=10) %>%
  mutate(component = fct_inorder(component))

terms_levels <- 
  plot_df %>% 
  group_by(terms) %>% 
  summarise(sum_value=sum(abs(value))) %>% 
  arrange(sum_value) %>% 
  pull(terms)

plot_df$terms <- factor(plot_df$terms, levels=terms_levels)

ggplot(plot_df, aes(value, terms, fill = abs(value))) +
  geom_col() +
  
  scale_fill_viridis_c(name="Abs value") +
  
  facet_wrap(~component, nrow = 1) +
  labs(title="Which variables contributed the most (together) to a particular PC?", 
       subtitle="This can be confusing. Make sure to ask your class teacher about any doubts you might have",
       captions="Figure 4. Contribution of most important\n features to each component",
       y = NULL) +
  # Prettify plot a bit
  theme_bw() +
  theme(plot.title=element_text(size=rel(1.2)),
        plot.subtitle = element_text(size=rel(1)),
        axis.title=element_text(size=rel(1.3)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        axis.text.y=element_text(size=rel(1.25)))
```

🗣️ **Discussion**:

- How does Figure 4 help you interpret Figures 2 & 3?
- How does the above help you think about the attributes of the **most common type** of OK Cupid profiles?
- How would you interpret the 2D coordinates if instead of PC01 vs PC02, we had plotted PC01 vs PC03 in Figures 2 & 3?

> Add your answer and reflections here

## Part III: Anomaly detection techniques (40 mins)

We will try two different techniques for outlier detection on the two principal components: DBSCAN and LOF.

👥 IN PAIRS, go through the action points below and discuss your impressions and key takeaways.

🎯 **ACTION POINTS**

1. Take a look at the clusters identified by DBSCAN:

```{r fig.width=10, out.width=8, echo=FALSE}
# DBSCAN
eps = 1
minPts = 20

dbscan_fit <- 
  pca_results %>%
  select(PC01, PC02) %>%
  dbscan(eps=eps, minPts = minPts)

# Save the cluster back to the data frame
pca_results['dbscan_cluster'] = factor(dbscan_fit$cluster)

ggplot(pca_results, 
       aes(PC01, PC02, color=dbscan_cluster, alpha=dbscan_cluster)) +
  geom_point(size=rel(3.5)) +
  
  scale_colour_manual(name="DBSCAN\nCluster", values=c(`0`="#C63C4A", `1`="#878787")) +
  scale_alpha_manual(name="DBSCAN\nCluster", values=c(`0`=0.6, `1`=0.15)) +
  scale_y_continuous(breaks=seq(-7.5, 7.5, 2.5), limits=c(-7.5, 7.5)) +
  
  labs(title="DBSCAN helps us identify outliers", 
       subtitle="The algorithm identified 2 different clusters, labelled 0 & 1",
       captions="Figure 5. Clusters identified by DBSCAN") +
  # Prettify plot a bit
  theme_bw() +
  theme(plot.title=element_text(size=rel(1.5)),
        plot.subtitle = element_text(size=rel(1.2)),
        axis.title=element_text(size=rel(1.3)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        axis.text=element_text(size=rel(1.25)))
```

🗣 **Discussion**: How well do you think DBSCAN performs at anomaly detection on the two principle components?

> Write your answer/reflections here.

2. Let's filter to take a look at these samples:

(Ask your class teacher for help if you don't know what we are doing here)

```{r}
df_anomalies <-
  pca_results %>%
  mutate(outlier_type=case_when(
    dbscan_cluster == 0 & PC02 > 2.5 ~ "upper side",
    dbscan_cluster == 0 & PC02 < -2.5 ~ "lower side",
    .default="Not an outlier"
  )) %>%
  filter(outlier_type != "Not an outlier") %>%
  arrange(desc(outlier_type)) %>%
  select(outlier_type, PC01, PC02, age, height, sex_m)

df_anomalies
```

The problem with the data frame above is that we are looking at the standardised values. It will probably be more helpful if we look at the original data instead.

💡 Do you see why we selected these particular variables here?

3. **Reverse the normalisation:**

What were the mean and standard deviation of the normalized variables?

```{r}
normalisation_step <- pca_recipe %>% tidy() %>% filter(type == "normalize") %>% pull(number)
normalisation_stats <- tidy(pca_recipe, normalisation_step)
normalisation_stats
```
Or in a more useful format:

```r
normalisation_stats <- normalisation_stats %>% 
  pivot_wider(names_from="statistic", values_from="value")

normalisation_stats
```

Reverse-engineering the process for these variables:

```{r, echo=FALSE}
age_stats <- normalisation_stats %>% filter(terms == "age")
height_stats <- normalisation_stats %>% filter(terms == "height")

df_anomalies$age <- (df_anomalies$age * age_stats$sd) + age_stats$mean
df_anomalies$height <- (df_anomalies$height * height_stats$sd) + height_stats$mean

# Height is shown in inches, to see it in cm do df_anomalies %>% mutate(height=2.54*height)
df_anomalies
```

💡 These are unusually tall and unusually short people on the upper and lower sides of the outliers, respectively. 

4. Small digression off-topic: Compare it to how k-means divides the data into k=2 clusters:

```{r  fig.width=10, out.width=8, echo=FALSE}
kmeans_fit <-
  pca_results %>% 
  select(PC01, PC02) %>% 
  kmeans(centers=2)

# Save the cluster back to the data frame
pca_results['kmeans_cluster'] = factor(kmeans_fit$cluster)

ggplot(pca_results, 
       aes(PC01, PC02, color=kmeans_cluster)) +
  geom_point(size=rel(3.5), alpha=0.3) +
  
  scale_colour_discrete(name="k-means\nCluster") + 
  scale_y_continuous(breaks=seq(-7.5, 7.5, 2.5), limits=c(-7.5, 7.5)) +
  
  labs(title="k-means captures a completely different pattern", 
       subtitle="What do you make of it?",
       captions="Figure 6. Clusters identified by\nk-means with k=2 clusters") +
  # Prettify plot a bit
  theme_bw() +
  theme(plot.title=element_text(size=rel(1.5)),
        plot.subtitle = element_text(size=rel(1.2)),
        axis.title=element_text(size=rel(1.3)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        axis.text=element_text(size=rel(1.25)))

```

5. Let's see if the Local Outlier Factor (LOF) performs better than DBSCAN. 

We use the `lof()` function to calculate local outlier factors:  

```{r}
lof_results <- pca_results %>%
  select(PC01, PC02) %>%
  # 5 is number of neighbors used to calculate the LOF
  lof(5)

# Append the LOF score as a new column
pca_results$score <- lof_results
```

What is the distribution of scores?

```{r  out.width=7, out.width=8, echo=FALSE}
ggplot(pca_results, aes(x=score, fill=score >= 1.5)) +
  geom_histogram(binwidth=0.025) +
  
  scale_fill_manual(name="Comparatively\nHigh scores", 
                    values=c(`FALSE`="#878787", `TRUE`="#C63C4A")) +
  
  labs(title="The majority of samples have a score around 1.0", 
       subtitle="As we are looking for anomalies, our interest lies in the samples with outlier scores",
       captions="Figure 7. Distribution of LOF scores") +
  # Prettify plot a bit
  theme_bw() +
  theme(plot.title=element_text(size=rel(1.5)),
        plot.subtitle = element_text(size=rel(1.2)),
        axis.title=element_text(size=rel(1.3)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        axis.text=element_text(size=rel(1.25)),
        legend.position = c(0.85, 0.3))
```

We can now plot the results of LOF.

```{r  fig.width=10, out.width=8, echo=FALSE}
ggplot(pca_results, aes(PC01, PC02, size=score, alpha=score >= 1.5, color=score)) +
  geom_point() +

  scale_alpha_manual(name="LOF scores", values=c(`FALSE`=0.1, `TRUE`=0.6), guide="none") +
  scale_radius(name="LOF scores", range=c(1, 10)) +
  scale_color_viridis_c(name="LOF scores") +  
  
    labs(title = "Samples coloured and sized according to their LOF scores", 
       caption = "Figure 8. Samples coloured and sized\naccording to their LOF scores") +
  

  # I tweaked those manually after seeing the plot
  scale_y_continuous(breaks=seq(-7.5, 7.5, 2.5), limits=c(-7.5, 7.5)) +
  
  # Prettify plot a bit
  theme_bw() +
  theme(plot.title=element_text(size=rel(1.2)),
        plot.subtitle = element_text(size=rel(1)),
        axis.title=element_text(size=rel(1.3)),
        axis.title.x=element_text(margin=margin(t=10)),
        axis.title.y=element_text(margin=margin(r=10)),
        axis.text=element_text(size=rel(1.2))) 
```

🗣 **Discussion**: Does LOF perform better than DBSCAN to detect 'anomalous' samples?  


GENIAL Participants:

-   **Fill out this brief survey at the end of the lab:** 🔗 [link](https://forms.office.com/e/h0dXriciyy) (requires LSE login)
  - You don't have to use ChatGPT today, but if you do, please share the chat log!