---
title: "Code for Session 8"
author: "Dr. Richard M. Crowley"
format:
  html:
    code-tools: true
    code-link: true
    embed-resources: true
---

Note that the directories used to store data are likely different on your computer, and such references will need to be changed before using any such code.

```{r helpers, warning=FALSE, message=F}
library(knitr)
library(kableExtra)

#' Pretty formatted HTML tables
#' 
#' @param df The dataframe to format as a table
#' @param names Names to use for columns if not using the df column names. Vector
#' @param highlight_cols Columns to highlight. Default is c(1) for highlighting the first row. Use c() to remove.
#' @param highlight_rows Rows to highlight. Default is c() to not highlight any rows.
#' @param color_rows Rows to highlight in color. Default is c() to not highlight any rows.
#' @param highlight_end_rows Rows to highlight. Uses a reverse indexing where 1 is the last row, 2 is second to last, etc.
#' @param full Use to span the full width of a slide. FALSE by default.
#' @param fixed_header Whether the first row should be treated as a header row. TRUE by default.
html_df <- function(df, names=NULL, highlight_cols=c(1), highlight_rows=c(), color_rows=c(), highlight_end_rows=c(), full=F, fixed_header=T) {
  if(!length(names)) {
    names=colnames(df)
  }
  kbl(df,"html", col.names = names, align = c("l",rep('c',length(cols)-1))) %>%
    kable_paper(c("striped","hover"), full_width=full, fixed_thead=F, html_font = "\"Source Sans Pro\", Helvetica, sans-serif") %>%
    {if(length(highlight_cols)) column_spec(., highlight_cols,bold=T, background = "#ffffff99") else .} %>%
    {if(length(highlight_rows)) row_spec(., highlight_rows,bold=T, background = "#ffffff99") else .} %>%
    {if(length(color_rows)) row_spec(., color_rows,bold=T, background = "#aaaaff99") else .} %>%
    {if(length(highlight_end_rows)) row_spec(., nrow(df) + 1 - highlight_end_rows,bold=T, background = "#ffffff99") else .}
}

library(tidyverse)
```

```{r}
#| include: false
options(downlit.attached = c("knitr", "kableExtra", "tidyverse"))
```

```{r}
readRDS("../../Data/Session_8-corp_summary.rds")
```

```{r}
corp_FOG <- readRDS("../../Data/Session_8-corp_readability_FOG.rds")
corp_FOG %>%
  head() %>%
  html_df()
```

```{r, fig.height=4}
summary(corp_FOG$FOG)
ggplot(corp_FOG, aes(x=FOG)) + geom_density() + xlim(c(10,40))
```

```{r, warning=F}
df_SIC <- read.csv('../../Data/Session_8-Filings2021.csv') %>%
  select(accession, regsic) %>%
  mutate(accession=paste0(accession, ".txt")) %>%
  rename(document=accession) %>%
  mutate(industry = case_when(
    regsic >=0100 & regsic <= 0999 ~ "Agriculture",
    regsic >=1000 & regsic <= 1499 ~ "Mining",
    regsic >=1500 & regsic <= 1799 ~ "Construction",
    regsic >=2000 & regsic <= 3999 ~ "Manufacturing",
    regsic >=4000 & regsic <= 4999 ~ "Utilities",
    regsic >=5000 & regsic <= 5199 ~ "Wholesale Trade",
    regsic >=5200 & regsic <= 5999 ~ "Retail Trade",
    regsic >=6000 & regsic <= 6799 ~ "Finance",
    regsic >=7000 & regsic <= 8999 ~ "Services",
    regsic >=9100 & regsic <= 9999 ~ "Public Admin" )) %>%
  group_by(document) %>%
  slice(1) %>%
  ungroup()
```

```{r, warning=F}
corp_FOG <- corp_FOG %>% left_join(df_SIC)
```

```{r}
corp_FOG %>%
  head() %>%
  html_df()
```

```{r, fig.height=5, warning=F, message=F}
ggplot(corp_FOG[!is.na(corp_FOG$industry),], aes(x=factor(industry), y=FOG)) +
  geom_violin(draw_quantiles = c(0.25, 0.5, 0.75)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + ylim(c(10, 40))
```

```{r, fig.height=4}
ggplot(corp_FOG[!is.na(corp_FOG$industry),], aes(x=FOG)) +
  geom_density() + facet_wrap(~industry) + xlim(c(10, 40))
```

```{r, message=F, warning=F}
library(DT)
readRDS('../../Data/Session_8-corp_kwic_gw.rds') %>%
  as.tibble() %>%
  mutate(text=paste(pre,keyword,post)) %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  select(docname, industry, text) %>%
  sample_n(100) %>%
  datatable(options = list(pageLength = 5), rownames=F)
```

```{r, message=F, warning=F}
gw_count <- readRDS('../../Data/Session_8-corp_kwic_gw.rds') %>%
  as.tibble() %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  group_by(docname) %>%
  slice(1) %>%
  ungroup() %>%
  group_by(industry) %>%
  mutate(docs_gw = n()) %>%
  slice(1) %>%
  ungroup() %>%
  select(industry, docs_gw)

corp_FOG %>% group_by(industry) %>%
  mutate(docs = n()) %>%
  slice(1) %>%
  ungroup() %>%
  left_join(gw_count) %>%
  mutate(docs_gw = ifelse(is.na(docs_gw),0,docs_gw)) %>%
  mutate(`Global warming` = docs_gw / docs,
         `Does not mention` = (docs - docs_gw) / docs) %>%
  gather(mention, percent, `Global warming`, `Does not mention`) %>%
  ggplot() + 
  geom_col(aes(x = industry, y = percent, fill=mention)) + 
  ylab("Percent of annual reports") +
  xlab("Industry (SIC code)") +
  ggtitle("Industries discussing global warming in 2021") + 
  scale_y_continuous(labels = scales::percent) + 
  scale_fill_manual(values=c("#CCCCCCAA", "#88CC88")) +
  theme(axis.text.x = element_text(angle=45, hjust=1))
```

```{r, message=F, warning=F}
library(DT)
readRDS('../../Data/Session_8-corp_kwic_covid.rds') %>%
  as.tibble() %>%
  mutate(text=paste(pre,keyword,post)) %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  select(docname, industry, text) %>%
  sample_n(100) %>%
  datatable(options = list(pageLength = 5), rownames=F)
```

```{r, message=F, warning=F}
covid_count <- readRDS('../../Data/Session_8-corp_kwic_covid.rds') %>%
  as.tibble() %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  group_by(docname) %>%
  slice(1) %>%
  ungroup() %>%
  group_by(industry) %>%
  mutate(docs_covid = n()) %>%
  slice(1) %>%
  ungroup() %>%
  select(industry, docs_covid)

corp_FOG %>% group_by(industry) %>%
  mutate(docs = n()) %>%
  slice(1) %>%
  ungroup() %>%
  left_join(covid_count) %>%
  mutate(docs_covid = ifelse(is.na(docs_covid),0,docs_covid)) %>%
  mutate(`COVID` = docs_covid / docs,
         `Does not mention` = (docs - docs_covid) / docs) %>%
  gather(mention, percent, `COVID`, `Does not mention`) %>%
  ggplot() + 
  geom_col(aes(x = industry, y = percent, fill=mention)) + 
  ylab("Percent of annual reports") +
  xlab("Industry (SIC code)") +
  ggtitle("Industries discussing COVID in 2021") + 
  scale_y_continuous(labels = scales::percent) + 
  scale_fill_manual(values=c("#CCCCCCAA", "#88CC88")) +
  theme(axis.text.x = element_text(angle=45, hjust=1))
```

```{r}
readRDS('../../Data/Session_8-corp_dfm_feat.rds')
```

```{r}
readRDS('../../Data/Session_8-corp_tfidf_feat.rds')
```

```{r}
out <- readRDS('../../Data/Session_8-corp_out_stm.rds')
#out <- readRDS('../../Data/Session_8-corp_out_lda.rds')
#out <- readRDS('../../Data/Session_8-corp_out_topicmodels.rds')
```

```{r}
out$documents[[1]][,386:390]
out$vocab[c(out$documents[[1]][,386:390][1,])]
```

```{r}
library(stm)
topics <- readRDS('../../Data/Session_8-corp_stm_topics.rds')
```

```{r}
labelTopics(topics)
```

```{r, message=F}
out$meta$industry <- factor(out$meta$industry)

doc_topics = data.frame(document=names(out$documents),
                        industry=out$meta$industry,
                        topic=1,
                        weight=topics$theta[,1])
for (i in 2:10) {
  temp = data.frame(document=names(out$documents),
                    industry=out$meta$industry,
                    topic=i,
                    weight=topics$theta[,i])
  doc_topics = rbind(doc_topics, temp)
}
# Proportional topics (%)
doc_topics <- doc_topics %>%
  group_by(document) %>%
  mutate(topic_prop = weight / sum(weight)) %>%
  ungroup()
```

```{r, message=F}
# Manually label topics
topic_labels = data.frame(topic = 1:10,
  topic_name = c('XBRL', 'Banking', 'Services', 'Equity',
                 'Investment', 'Energy', 'R&D',
                 'Compensation', 'Financial', 'Debt'))

doc_topics <- doc_topics %>% left_join(topic_labels)
```

```{r}
doc_topics %>% filter(document=='0001564590-21-039151.txt')
```

```{r, fig.height=4}
doc_topics %>%
  filter(document=='0001564590-21-039151.txt' |
         document=='0000320193-21-000105.txt') %>%
  mutate(Company=ifelse(document=='0001564590-21-039151.txt', 'MSFT','AAPL')) %>%
  ggplot(aes(x=factor(topic_name), y=topic_prop, fill=factor(topic_name))) + 
  geom_col() + facet_wrap(~Company) + 
  theme(axis.text.x=element_blank(),axis.ticks.x = element_blank())
```

```{r, fig.height=4, warning=F}
doc_topics %>%
  group_by(industry, topic) %>%
  mutate(topic_prop = mean(topic_prop)) %>%
  slice(1) %>%
  ungroup() %>%
  ggplot(aes(x=factor(topic_name), y=topic_prop, fill=factor(topic_name))) + 
  geom_col() + facet_wrap(~industry) +
  theme(axis.text.x=element_blank(),axis.ticks.x = element_blank())
```

```{r, message=F}
library(tidyr)
wide_topics <- spread(doc_topics[,c(1,2,5,6)], topic_name, topic_prop)
# Note: dropping XBRL here
mat <- wide_topics[,3:11]

mat[,1:6] %>% head(n=3) %>% html_df(highlight_cols = c())
```

```{r}
set.seed(6845868)
clusters <- kmeans(mat, 9)

# Add clusters back into our data
wide_topics$kmean <- clusters$cluster

clusters$cluster %>% head()
```

```{r, fig.height=4}
cbind(as.data.frame(clusters$center), data.frame(kmean=1:9)) %>%
  gather("Topics","weights",-kmean) %>%
  ggplot(aes(x=factor(Topics), y=weights, fill=factor(Topics))) +
  geom_col() + 
  facet_wrap(~kmean) + 
  theme(axis.text.x=element_blank(),axis.ticks.x = element_blank())
```

```{r}
#| eval: true
#| echo: true

library(uwot)

# Build the UMAP model
umap_train <- umap(mat, ret_model = TRUE)
# Extract coordinates
umap_coords <- umap_train$embedding %>% as.data.frame()
colnames(umap_coords) <- c('umap1', 'umap2')
# Merge coordinates into our data frame
wide_topics <- cbind(wide_topics, umap_coords)
```

```{r, fig.height=5}
ggplot(wide_topics,
       aes(x = umap1, y = umap2,
           color = industry)) + 
  geom_point(alpha = 0.3) +
  theme_bw()
```

```{r, fig.height=5}
ggplot(wide_topics,
       aes(x = umap1, y = umap2,
           color = factor(kmean))) + 
  geom_point(alpha = 0.3) +
  theme_bw()
```

```{r, fig.height=5}
ggplot(wide_topics, aes(x=kmean)) + geom_bar() + facet_wrap(~factor(industry))
```

```{r, fig.height=5}
ggplot(wide_topics, aes(x=umap1, y=umap2, color=factor(kmean))) + geom_point() +
  facet_wrap(~factor(industry))
```

```{r, fig.height=5}
ggplot(wide_topics, aes(x=umap1, y=umap2, color=factor(industry))) + geom_point() +
  facet_wrap(~factor(kmean))
```

```{r}
wide_topics$dist <- sqrt(rowSums(abs(mat - fitted(clusters))))  # Distance from center
wide_topics[,c(1,2,4,8,16)] %>% arrange(desc(dist)) %>% slice(1:5) %>% html_df()
```

```{r}
wide_topics[,c(1,2,4,8,16)] %>%
  filter(industry!="Finance") %>%
  arrange(desc(dist)) %>%
  mutate(id=1:n()) %>%
  select(id,everything()) %>%
  slice(1:7) %>%
  html_df()
```

```{r}
wide_topics[,c(1,2,4,8,16)] %>% filter(Compensation > 0.9, dist > 1.1) %>% nrow()
```

```{r}
wide_topics[,c(1,2,4,8,16)] %>%
  filter(industry!="Finance", Compensation < 0.5) %>%
  arrange(desc(dist)) %>%
  mutate(id=1:n()) %>%
  select(id,everything()) %>%
  slice(1,2,3,8,9,10) %>%
  html_df()

```

```{r, message=F}
train <- wide_topics %>% filter(!is.na(industry))
label <- wide_topics %>% filter(is.na(industry))
```

```{r, message=F}
library(caret)
tout <- readRDS('../../Data/Session_8-corp_knn.rds')
```

```{r}
tout
```

```{r, fig.height=4}
ggplot(tout$results, aes(x=k, y=Accuracy)) +
  geom_line() + 
  geom_ribbon(aes(ymin=Accuracy - AccuracySD*1.96,
                  ymax=Accuracy + AccuracySD*1.96), alpha=0.2) + 
  geom_vline(xintercept=15, color="blue") + 
  xlab("k, optimal = 15")
```

```{r}
label$industry_pred <- predict(tout,
                              label)
label[,c("document",
         "industry_pred")] %>%
  head %>% html_df
```

