---
title: "Code for Session 8"
author: "Dr. Richard M. Crowley"
date: ""
output:
  html_notebook
---

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}
library(knitr)
library(kableExtra)
html_df <- function(text, cols=NULL, col1=FALSE, full=F) {
  if(!length(cols)) {
    cols=colnames(text)
  }
  if(!col1) {
    kable(text,"html", col.names = cols, align = c("l",rep('c',length(cols)-1))) %>%
      kable_styling(bootstrap_options = c("striped","hover"), full_width=full)
  } else {
    kable(text,"html", col.names = cols, align = c("l",rep('c',length(cols)-1))) %>%
      kable_styling(bootstrap_options = c("striped","hover"), full_width=full) %>%
      column_spec(1,bold=T)
  }
}
```

```{r}
library(tidyverse)
```

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

```{r}
corp_FOG <- readRDS("../../Data/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()
```

```{r, warning=F}
df_SIC <- read.csv('../../Data/Filings2014.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()
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))
```

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

```{r, fig.height=4}
library(lattice)
densityplot(~FOG | industry,
            data=corp_FOG,
            plot.points=F,
            main="Fog index distibution by industry (SIC)",
            xlab="Fog index",
            layout=c(3,3))
```

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

```{r, message=F, warning=F}
gw_count <- readRDS('../../Data/corp_kwic.rds') %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  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 2014") + 
  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/corp_dfm_feat.rds')
```

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

```{r}
readRDS('../../Data/corp_tfidf_bank.rds')
```

```{r}
out <- readRDS('../../Data/corp_out_stm.rds')
#out <- readRDS('../../Data/corp_out_lda.rds')
#out <- readRDS('../../Data/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/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('Real Estate', 'Management', 'Product',
                                         'Investment', 'Services', 'Financing',
                                         'Service2', 'Insurance', 'Industrial',
                                         'Utility'))

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

```{r}
doc_topics %>% filter(document=='0001104659-14-015152.txt')
```

```{r, fig.height=4}
doc_topics %>%
  filter(document=='0001104659-14-015152.txt' |
         document=='0000019617-14-000289.txt') %>%
  mutate(Company=ifelse(document=='0001104659-14-015152.txt', 'Citi','JPM')) %>%
  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)
mat <- wide_topics[,3:12]

mat[,1:6] %>% head() %>% html_df()
```

```{r}
set.seed(6845868)
clusters <- kmeans(mat, 9)
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, fig.height=4}
library(cluster) # Uses PCA (principle component analysis)
clusplot(mat, clusters$cluster, color=TRUE, shade=TRUE, 
         labels=4)
```

```{r, message=F, warning=F}
library(Rtsne)
dups <- which(duplicated(mat))
wide_nodup <- wide_topics[-dups,]
wide_nodup$kmean <- clusters$cluster[-dups]
tsne_data <- readRDS('../../Data/corp_tsne.rds')

wide_nodup <- wide_nodup %>%
  mutate(tsne1 = tsne_data$Y[, 1], tsne2 = tsne_data$Y[, 2])
```

```{r, fig.height=5}
ggplot(wide_nodup, aes(x = tsne1, y = tsne2, colour = industry)) + 
    geom_point(alpha = 0.3) + theme_bw()
```

```{r, fig.height=5}
ggplot(wide_nodup, aes(x = tsne1, y = tsne2, colour = factor(kmean))) + 
    geom_point(alpha = 0.3) + theme_bw()
```

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

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

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

```{r}
#wide_topics$dist <- sqrt(rowSums(mat - fitted(clusters))^2)
wide_topics$dist <- sqrt(rowSums(abs(mat - fitted(clusters))))
wide_topics[,c(1,2,3,5,13)] %>% arrange(desc(dist)) %>% slice(1:5) %>% html_df()
```

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

```{r}
wide_topics[,c(2,6,9,10,11,12,13)] %>% filter(industry!="Finance") %>% arrange(desc(dist)) %>% mutate(id=1:n())%>% select(id,everything()) %>% slice(5) %>% html_df()
```

```{r, message=F}
train <- wide_topics[!is.na(wide_topics$industry),]
label <- wide_topics[is.na(wide_topics$industry),]
```

```{r, message=F}
library(caret)
tout <- readRDS('../../Data/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
```

```{r, message=F, warning=F}
ts_wt <- wide_nodup %>% left_join(label[,c("document","industry_pred")])

ts_wt <- ts_wt %>%
  mutate(tsne1 = tsne_data$Y[, 1], tsne2 = tsne_data$Y[, 2])

# Force consistent factor values
inds <- unique(ts_wt$industry)
ts_wt$industry <- factor(ts_wt$industry, inds)
ts_wt$industry_pred <- factor(ts_wt$industry_pred, inds)

# Replicate default ggplot colors
ggplotColours <- function(n = 6, h = c(0, 360) + 15){
  if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
  hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}

ggplot() +
  scale_shape_identity() + # Allow for more plot point options
  geom_point(data=ts_wt[!is.na(ts_wt$industry),],
             aes(x=tsne1, y=tsne2, color=industry, shape=1), size=1) + 
  geom_point(data=ts_wt[!is.na(ts_wt$industry_pred),], aes(x=tsne1, y=tsne2,
             fill=industry_pred, shape=23, stroke=0.5), size=2) +
  guides(fill = "none") + 
  scale_color_manual(values=ggplotColours(n = 9), labels=inds, drop=FALSE) + 
  scale_fill_manual(values=ggplotColours(n = 9), labels=inds, drop=FALSE)
```

```{r}
clusters <- kmeans(mat, 20)
clusters$cluster %>% head()
```

```{r, fig.height=4}
wide_nodup$kmean2 <- clusters$cluster[-dups]
ggplot(wide_nodup, aes(x = tsne1, y = tsne2, colour = factor(kmean2))) + 
    geom_point(alpha = 0.3) + theme_bw()
```

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

```{r, fig.height=5}
ggplot(wide_nodup, aes(x=tsne1, y=tsne2, color=factor(kmean2))) + geom_point() +
  facet_wrap(~factor(industry))
```

