---
title:
output:
  revealjs::revealjs_presentation:
    theme: white
    css: white-edits.css
    highlight: tango
    transition: slide
    self_contained: false
    reveal_plugins: ["chalkboard","search"]
    transitionSpead: fast
    mathjax: "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_SVG"
    reveal_options:
      minScale: 1.0
      maxScale: 1.0
      width: 1024
      height: 768
      slideNumber: true
      controlsTutorial: true
      chalkboard:
        theme: whiteboard
        readOnly: false
---

## ACCT 420: Topic modeling and anomaly detection {data-background="../Backgrounds/title.jpg" class="header no-toc-progress" data-state="no-toc-progress"}

<!-- "reveal_plugins:..." and "chalkboard:" are handled by prep.py -->
```{r helpers, echo=FALSE, 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)
  }
}
```


<h3 class="date"><br>Session 8</h3>
<h2 class="author"><br>Dr. Richard M. Crowley</h2>

<!-- Load in primary data set -->
```{r, include=FALSE}
library(tidyverse)
```

# Front matter {data-background="../Backgrounds/section.jpg" class="header"}

## Learning objectives {data-background="../Backgrounds/default.jpg" class="default"}

<div class="clearfix">
<div class="left">
<img src="../Course_Diagrams/Map-8.png">
</div>
<div class="right objectives">
- <b>Theory:</b>
    - NLP
    - Anomaly detection
- <b>Application:</b>
    - Understand annual report readability
    - Examine the *content* of annual reports
    - Group firms on content
    - Fill in missing data
- <b>Methodology:</b>
    - ML/AI (LDA)
    - ML/AI (k-means, t-SNE)
    - More ML/AI (KNN)
</div>
</div>

## Datacamp {data-background="../Backgrounds/default.jpg" class="default"}

> - One last chapter: <a href="https://www.datacamp.com/courses/introduction-to-machine-learning-with-r">What is Machine Learning</a>
>     - Just the first chapter is required
>     - You are welcome to do more, of course
> - This is the last *required* chapter on Datacamp

## Group project {data-background="../Backgrounds/presentation.jpg" class="default"}

> For reading large files, `package:readr` is your friend

```{r, eval=F}
library(readr)  # or library(tidyverse)
df <- read_csv("really_big_file.csv.zip")
```

- It can read directly from zip files!
    - Like those that you can export from WRDS
    - Good for saving disk space
- It can write directly to zip files too

## Group project {data-background="../Backgrounds/presentation.jpg" class="default"}


> For saving intermediary results, `saveRDS()` + `readRDS()` is your friend

```{r, eval=F}
saveRDS(really_big_object, "big_df.rds")

# Later on...
df <- readRDS("big_df.rds")
```

- You can neatly save processed data, finished models, and more
    - This is particularly helpful if you want to work on something later or distribute data or results to teammates

> If you look at the code file for this lesson, you'll see this used extensively

# Sets of documents (corpus) {data-background="../Backgrounds/section.jpg" class="header"}

## Importing sets of documents (corpus) {data-background="../Backgrounds/reference.jpg" class="default"}

- I will use the `package:readtext` package for this example
    - Importing all 6,000 annual reports from 2014
- Other options include using
    - `package:purrr` and `df_map()`
    - `package:tm` and `VCorpus()`
    - `package:textreadr` and `read_dir()`

```{r, eval=F}
library(readtext)
library(quanteda)
# Needs ~1.5GB
corp <- corpus(readtext("/media/Scratch/Data/Parser2/10-K/2014/*.txt"))
```

## Corpus summary {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, eval=F}
summary(corp)
```
```{r, echo=F}
readRDS("../../Data/corp_summary.rds")
```

## Running readability across the corpus {data-background="../Backgrounds/default.jpg" class="default"}

```{r, eval=F}
# Uses ~20GB of RAM...  Break corp into chunks if RAM constrained
corp_FOG <- textstat_readability(corp, "FOG")
corp_FOG %>%
  head() %>%
  html_df()
```

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

> Recall that Citi's annual report had a Fog index of 21.63

## Readability across documents {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

## Are certain industries' filings more readable? {data-background="../Backgrounds/default.jpg" class="default"}

- Since the SEC has their own industry code (SIC), we'll use <a href="https://www.osha.gov/pls/imis/sic_manual.html">that</a>
- SIC codes are 4 digits
    - The first two digits represent the industry
    - The third digit represents the business group
    - The fourth digit represents the specialization
- Example: Citigroup is SIC 6021
    - 60: Depository institution
    - 602: Commercial bank
    - 6021: National commercial bank

## Are certain industries' filings more readable? {data-background="../Backgrounds/default.jpg" class="default"}

- Merge in SIC code by group

```{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)
```

## Are certain industries' filings more readable? {data-background="../Backgrounds/default.jpg" class="default"}

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

## Are certain industries' filings more readable? {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{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))
```

## Are certain industries' filings more readable? {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

## Are certain industries' filings more readable? {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{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))
```

## `package:quanteda` bonus: Finding references across text {data-background="../Backgrounds/default.jpg" class="default"}

```{r, eval=F}
kwic(corp, phrase("global warming")) %>% mutate(text=paste(pre,keyword,post)) %>%
  select(docname, text) %>% datatable(options = list(pageLength = 5), rownames=F)
```
```{r, echo=F, 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)
```

## `package:quanteda` bonus: Mentions by industry {data-background="../Backgrounds/default.jpg" class="default"}

```{r, echo=F, 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))
```

# Going beyond simple text measures {data-background="../Backgrounds/section.jpg" class="header"}

## What's next {data-background="../Backgrounds/default.jpg" class="default"}

- Armed with an understanding of how to process unstructured data, all of the sudden the amount of data available to us is expanding rapidly
- To an extent, anything in the world can be viewed as data, which can get overwhelming pretty fast
- We'll require some better and newer tools to deal with this

## Problem: What do firms discuss in annual reports? {data-background="../Backgrounds/case.jpg" class="default"}

- This is a hard question to answer -- our sample has 104,690,796 words in it!
    - 69.8 hours for the "world's fastest reader", per <a href="https://skeptoid.com/episodes/4229">this source</a>
    - 103.86 **days** for a standard speed reader (<a href="http://psycnet.apa.org/record/1998-11174-004">700wpm</a>)
    - 290.8 days for an average reader (<a href="http://psycnet.apa.org/record/1998-11174-004">250wpm</a>)
- We could read a small sample of them?
- Or... have a computer read all of them!

## Recall the topic variable from session 6 {data-background="../Backgrounds/default.jpg" class="default"}

- Topic was a set of 31 variables indicating *how much* a given topic was discussed
- This measure was created by making a machine read every annual report
    - The computer then used a technique called LDA to process these reports' content into topics

<img src="../Figures/artificial-neural-network-3501528_640.png" height="200px">

> This is our end goal, but we'll work our way up

## Term document matrices (TDM) {data-background="../Backgrounds/default.jpg" class="default"}

- Before we begin, we'll need a matrix of word counts per document
- We'll create something called a *sparse matrix* for this
- A sparse matrix is a matrix that only lists values that aren't 0

> Think about the structure of a matrix where rows are document names and columns are individual words.  How much of this matrix will be 0s?

## Making a TDM {data-background="../Backgrounds/reference.jpg" class="default"}

- In `package:quanteda`, use `dfm()`
    - Useful options:
        - `stem=TRUE`, Code similar words as the same
            - Ex.: *cod*e, *cod*ing, and *cod*er would all become *cod*
                - Helps with the curse of dimensionality
        - `remove=c(...)`, You can supply a list of stop words to remove
            - You can use `remove=stopwords()` for a simple list
            - The `stopwords()` function is provided by the `package:stopwords` package, and actually supports over 50 languages, including Chinese, English, Hindi, and Malay
            - We can use SMART like last week: `remove=stopwords(source='smart')`
            - For other languages, use `remove=stopwords("zh", source="stopwords-iso")`

## Making a TDM {data-background="../Backgrounds/default.jpg" class="default"}

```{r, eval=F}
# adding industry to the corpus
docs <- docnames(corp)
docs <- data.frame(document=docs)
docs <- docs %>% left_join(df_SIC)
docvars(corp, field="industry") <- docs$industry
```

```{r, eval=F}
# Simplest way
tdm <- dfm(corp)

# With stopwords
tdm <- dfm(corp,
           remove=stopwords(source='smart'))

# With stopwords and stemming -> Used in next slides
# 2.6B elements in the matrix
tdm <- dfm(corp,
           stem=TRUE,
           remove=stopwords(source='smart'),
           remove_punct=TRUE,
           remove_numbers=TRUE) %>%
  dfm_trim(min_termfreq=10, termfreq_type = "count")
```

## What words matter by industry? {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, eval=F}
topfeatures(tdm, n=5, groups="industry")
```
```{r, echo=F}
readRDS('../../Data/corp_dfm_feat.rds')
```

> This isn't very informative

## TF-IDF {data-background="../Backgrounds/reference.jpg" class="default"}

- Words counts are not very informative
- Knowing the words that show up frequently in one group *but not in the others* would be much more useful
- This is called TF-IDF
    - **T**erm **F**requency-**I**nverse **D**ocument **F**requency
- Think of it roughly as:

$$
\frac{\text{How many times a word is in the document}}{\text{How many documents the word is in}}
$$

- We can easily calculate TF-IDF using `dfm_tfidf()` from `package:quanteda`
    - The options we'll specify are used to match a more standard output
    - `package:quanteda`'s default options are a bit odd

## The actual TF-IDF equation we'll use {data-background="../Backgrounds/default.jpg" class="default"}

$$
 \frac{f_{w,d}}{f_d} \cdot -\log_2\left(\frac{n_t}{N}\right)
$$

- $w$ represents 1 word
- $d$ represents 1 document
- $f_{w,d}$ is the number of times $w$ appears in $d$
- $f_{d}$ is the number of times any word appears in $d$
- $n_w$ is the number of documents with $w$ at least once
- $N$ is the number of documents


## What words matter by industry? {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, eval=F}
tfidf_mat <- dfm_tfidf(tdm, base=2, scheme_tf="prop")
topfeatures(tfidf_mat, n=5, groups="industry")
```

```{r, echo=F}
readRDS('../../Data/corp_tfidf_feat.rds')
```

> These are more meaningful

## What words matter for banks? {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r,eval=F}
topfeatures(tfidf_mat, n=20, groups="industry")$Finance
```
```{r, echo=F}
readRDS('../../Data/corp_tfidf_bank.rds')
```

# Moving on to LDA {data-background="../Backgrounds/section.jpg" class="header"}

## What is LDA? {data-background="../Backgrounds/reference.jpg" class="default"}

- **L**atent **D**irichlet **A**llocation
- One of the most popular methods under the field of *topic modeling*
- LDA is a Bayesian method of assessing the content of a document
- LDA assumes there are a set of topics in each document, and that this set follows a *Dirichlet* prior for each document
    - Words within topics also have a *Dirichlet* prior

<img src="../Figures/dirichlet.png" height="100px">

> <a href="http://www.cs.columbia.edu/~blei/papers/Blei2012.pdf">More details from the creator</a>

## An example of LDA {data-background="../Backgrounds/default.jpg" class="default"}

<img src="../Figures/LDA_blei.png" height="600px">

## How does it work? {data-background="../Backgrounds/default.jpg" class="default"}

1. Reads all the documents
    - Calculates counts of each word within the document, tied to a specific ID used across all documents
2. Uses variation in words within and across documents to infer topics
    - By using a Gibbs sampler to simulate the underlying distributions
        - An MCMC method

- It's quite complicated in the background, but it boils down to a system where generating a document follows a couple rules:
   1. Topics in a document follow a multinomial/categorical distribution
   2. Words in a topic follow a multinomial/categorical distribution
   
## Implementations in R {data-background="../Backgrounds/reference.jpg" class="default"}

- There are at least four good implementations of LDA in R
    1. `package:stm`: A bit of a tweak on the usual LDA model that plays nicely with `package:quanteda` and also has an associated `package:stmBrowser` package for visualization (on Github)
    2. `package:lda`: A somewhat rigid package with difficult setup syntax, but it plays nicely with the great `package:LDAvis` package for visualizing models.  Supported by `package:quanteda`.
    3. `package:topicmodels`: An extensible topic modeling framework that plays nicely with `package:quanteda`
    4. `package:mallet`: An R package to interface with the venerable <a href="http://mallet.cs.umass.edu/">MALLET Java package</a>, capable of more advanced topic modeling



## Implementing a topic model in STM {data-background="../Backgrounds/default.jpg" class="default"}

```{r, eval=F}
# quanteda's conversion for the stm package
out <- convert(tdm, to = 'stm')
# quanteda's conversion for the lda package
# out <- convert(tdm, to = 'lda')
# quanteda's conversion for the topicmodels package
# out <- convert(tdm, to = 'topicmodels')
```

```{r, echo=F}
out <- readRDS('../../Data/corp_out_stm.rds')
#out <- readRDS('../../Data/corp_out_lda.rds')
#out <- readRDS('../../Data/corp_out_topicmodels.rds')
```

- Creates a list of 3 items:
    - `out$documents`: Index number for each word with count/document
    - `out$vocab`: Words and their index numbers
    - `out$meta` a data frame of information from the corpus (`industry`)

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

## Running the model {data-background="../Backgrounds/default.jpg" class="default"}

- We will use the `stm()` function from the `package:stm` package
    - It has a lot of options that you can explore to tweak the model
    - The most important is `K`, the number of topics we want.  I'll use 10 for simplicity, but often we need more to neatly categorize the text
        - `K=100` is a popular choice when we are using the output as an input to another model
        - The model we used in Session 6 had `K=31`, as that captures the most restatements in sample

```{r, eval=F}
library(stm)
topics <- stm(out$documents, out$vocab, K=10)
```
```{r, include=F}
library(stm)
topics <- readRDS('../../Data/corp_stm_topics.rds')
```

> <a href="./corp_stm_print.txt">What this looks like while running</a>

## LDA model {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

- `Highest prob` is a straightforward measure to interpret
    - The words with the highest probability of being chosen in the topic

## Applying our topic model to our data {data-background="../Backgrounds/default.jpg" class="default"}

```{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)
```

## Topic content of the Citi 10-K {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

## Topic content of the Citi 10-K versus JPMorgan {data-background="../Backgrounds/default.jpg" class="default"}

```{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())
```

## Topic content by industry {data-background="../Backgrounds/default.jpg" class="default"}

```{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())
```

## A nice visualization of our STM model {data-background="../Backgrounds/microscope.jpg" class="default"}

- Using LDAvis via `package:STM`'s `toLDAvis()` function
    - Need `package:LDAvis` and `package:servr` installed to run

```{r, eval=F}
# Code to generate LDAvis
toLDAvis(topics, out$documents, R=10)
```

> <a href="LDAvis/index.html">Click to view</a>

- Using `package:stmBrowser`'s `stmBrowser()` function
    - Install from github

```{r, eval=F}
# code to generate stmBrowser
stmBrowser(topics, data=data.frame(text=names(out$documents),
                                   industry=out$meta$industry),
           c('industry'), text='text')
```

> <a href="stm-visualization/data/stm-visualization/index.html">Click to view</a>

## What we have accomplished? {data-background="../Backgrounds/default.jpg" class="default"}

- We have created a measure of the content of annual reports
    - This gives us some insight as to what is discussed in *any* annual report from 2014 by looking at only 10 numbers as opposed to having to read the whole document
        - We can apply it to other years as well, though it will be a bit less accurate if new content is discussed in those years
    - We can use this measure in a variety of ways
        - Some forecasting related, such as building in firm disclosure into prediction models
        - Some forensics related, such as our model in Session 6

## Consider {data-background="../Backgrounds/group.jpg" class="default"}

> How might we leverage LDA (or other topic modeling methods) to improve and simplify analytics?

- What other contexts or data could we use LDA on?
- What other problems can we solve with LDA?

# Clustering without known groups {data-background="../Backgrounds/section.jpg" class="header"}

## Problem: Classifying companies based on disclosure {data-background="../Backgrounds/case.jpg" class="default"}

- While industry code is one classification of firms, it has a number of drawbacks:
    1. The classification system is old and perhaps misses new industries
    2. It relies on self-reporting
    3. Firms' classifications rarely change, even when firms themselves change

> We'll build a different classification system, based on what they discuss in their annual reports

## Clustering {data-background="../Backgrounds/default.jpg" class="default"}

- One important aspect of detecting anomalies is determining groups in the data
    - We call this *clustering*
- If we find that a few elements of our data don't match the usual groups in the data, we can consider this to be an anomaly
    - Similar to the concept of outliers, but taking into account *multiple variables* simultaneously

<div class="clearfix">
<div class="left">
- The grey dot is at the mean of both the $x$ and $y$ dimensions
    - it isn't an outlier
- But there are 4 clear clusters... and it doesn't belong to any!
</div>
<div class="right">
<img src="../Figures/grouping_anom.png" height="250px">
</div>
</div>

## One clustering approach: k-means {data-background="../Backgrounds/default.jpg" class="default"}

$$
\min_{C_k} \sum_{k=1}^{K} \sum_{x_i\in C_k} \left(x_i - \mu_k\right)^2
$$

- Minimizes the sum of squared distance between points within groups
- Technically this is a machine learning algorithm, despite its simplicity
- You need to specify the number of groups you want

<div class="clearfix">
<div class="left">
- Pros:
    - Very fast to run
    - Simple interpretation
</div>
<div class="right">
- Cons
    - Simple algorithm
    - Need to specify $k$, the number of clusters
</div>
</div>

## Prepping data {data-background="../Backgrounds/default.jpg" class="default"}

- We will need data to be in a matrix format, with...
    - 1 row for each observation
    - 1 column for each variable we want to cluster by
- Since our data is currently in a long format, we'll recast this with `package:tidyr`


```{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()
```

## Calculating k-means {data-background="../Backgrounds/reference.jpg" class="default"}

```{r}
set.seed(6845868)
clusters <- kmeans(mat, 9)
clusters$cluster %>% head()
```

- The algorithm tells us group numbers for each observation
- The numbers themselves are arbitrary
    - The clustering (observations sharing a group number) is what matters

## Visualizing the clusters {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{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())
```

## Visualizing k-means with PCA {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, fig.height=4}
library(cluster) # Uses PCA (principle component analysis)
clusplot(mat, clusters$cluster, color=TRUE, shade=TRUE, 
         labels=4)
```

## Improving our visualization {data-background="../Backgrounds/default.jpg" class="default"}

- The PCA based map is really unreadable
    - This is usually the case, unless you have only a few dimensions to the data
- There is a relatively new method (2008), t-SNE, that is significantly better
    - **t**-distributed **S**tochastic **N**eighbor **E**mbedding
    - A machine learning algorithm designed to explain machine learning algorithms
        - It maintains neighbor relationships while reducing dimensions
    - It takes a much longer time to run than PCA, however
    - Implemented efficiently in R in the `package:Rtsne` package

## Visualizing with t-SNE: Running t-SNE {data-background="../Backgrounds/default.jpg" class="default"}

```{r, eval=F, message=F}
library(Rtsne)
dups <- which(duplicated(mat))
wide_nodup <- wide_topics[-dups,]
wide_nodup$kmean <- clusters$cluster[-dups]

#This is slow (it's O(n log(n))).  Original model was O(n^2) though
tsne_data <- Rtsne(mat[-dups,])

wide_nodup <- wide_nodup %>%
  mutate(tsne1 = tsne_data$Y[, 1], tsne2 = tsne_data$Y[, 2])
```
```{r, echo=F, 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])
```

## Visualizing with t-SNE: Industries {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

## Visualizing with t-SNE: k-means {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

## Why are these graphs different? {data-background="../Backgrounds/default.jpg" class="default"}

- Possibly due to...
    - Data: 10-K disclosure content doesn't fully capture industry inclusion
    - LDA: The measure is noisy -- it needs more data
    - SIC code: The measure doesn't cleanly capture industry inclusion
        - Some firms are essentially misclassified
- Recall, SIC covers Agriculture, Forestry and Fishing; Mining; Construction; Manufacturing; Transportation, Communications, Electric, Gas, and Sanitary Services; Wholesale Trade; Retail Trade; Finance, Insurance, and Real Estate; Services; Public Administration

## How related are clusters and industries? {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

## How related are clusters and industries? {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

## How related are clusters and industries? {data-background="../Backgrounds/case_digital.jpg" class="default"}

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

## Great examples of t-SNE usage {data-background="../Backgrounds/lab.jpg" class="default"}


- <a href="http://colah.github.io/posts/2014-10-Visualizing-MNIST/">Visualizing handwritten numbers</a>
- <a href="https://colah.github.io/posts/2015-01-Visualizing-Representations/big_vis/wiki.html">Visualizing Wikipedia articles</a>
    - <a href="http://colah.github.io/posts/2015-01-Visualizing-Representations/">The full blog post</a>, which is a great read about visualizing high-dimensional data

<a href="http://jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf"><img src="../Figures/tsne-MNIST-small.png" height="350px"></a>

## Looking for anomalies {data-background="../Backgrounds/case_digital.jpg" class="default"}

- k-means minimizes the distance from a central point
- We can look for the firms that are farthest from said point!

```{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()
```

- 5 standard insurance companies
    - SIC codes lump banks and insurance together, but they are actually very different industries
    - E.g.: <a href="https://www.sec.gov/Archives/edgar/data/1171500/000117150014000007/ptp2013_10k.htm">Platinum Underwriters Holdings</a>

## Looking for anomalies {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, echo=F}
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()
```

- 1-4, 9-10: Healthcare services + real estate (1: <a href="https://www.sec.gov/Archives/edgar/data/860730/000119312514070000/d660799d10k.htm">HCS Holdings</a>)
- 7: Healthcare and insurance management (<a href="https://www.sec.gov/Archives/edgar/data/1303313/000119312514086648/d646542d10k.htm">Magellan Health Services</a>)
- 6 & 8: Healthcare management (<a href="https://www.sec.gov/Archives/edgar/data/1035688/000104746914001304/a2218369z10-k.htm">Select Medical Holdings</a> & <a href="https://www.sec.gov/Archives/edgar/data/353230/000035323014000012/ocr10-k2013.htm">Omnicare</a>)

```{r, echo=F}
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()
```

- 5: Partnership for TV/internet/telco in 2 Sourthern US rural areas
    - <a href="https://www.sec.gov/Archives/edgar/data/843368/000119312514077529/d650448d10k.htm">Northland Cable Properties Eight Ltd. Ptr.</a>


## What we have accomplished {data-background="../Backgrounds/default.jpg" class="default"}

- We have created a classification of firms into discrete groups based on their disclosure content of their 10-K filings
    - The classification accounts for how similar each firm's content is to other firms' content
- We have used this classification to identify 10 firms which have non-standard accounting disclosures for their SIC code classification

> Text based industry classification using 10-Ks has been shown to be quite viable, such as in work by <a href="http://hobergphillips.usc.edu/industryclass.htm">Hoberg and Phillips.</a>

## Consider {data-background="../Backgrounds/group.jpg" class="default"}

> What else could we use clustering to solve?

- Where in business would we like to group something, but we don't know the groups?

# Filling in missing data {data-background="../Backgrounds/section.jpg" class="header"}

## Problem: Missing data {data-background="../Backgrounds/case.jpg" class="default"}

- You may have noticed that some of the `industry` measure was `NA`
- What if we want to assign an industry to these firms based on the content of their 10-K filings?

## Using k-means {data-background="../Backgrounds/default.jpg" class="default"}

- One possible approach we could use is to fill based on the category assigned by k-means
- However, as we saw, k-means and SIC code don't line up perfectly...
    - So using this classification will definitely be noisy

## A better approach with KNN {data-background="../Backgrounds/default.jpg" class="default"}

- KNN, or **K**-**N**earest **N**eighbors is a *supervised* approach to clustering
- Since we already have industry classifications for most of our data, we can use that structure to inform our assignment of the missing industry codes
- The way the model uses the information is by letting the nearest labeled points "vote" on what the point should be
    - Points are defined by 10-K content in our case

## Implementing KNN in R {data-background="../Backgrounds/default.jpg" class="default"}

- We'll use the `package:caret` package for this, as it will allow us to use k-fold cross validation to select a model
    - The same technique we used for LASSO and xgboost

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

```{r, eval=F}
library(caret)
trControl <- trainControl(method='cv', number=20)
tout <- train(industry ~ .,
      method = 'knn',
      tuneGrid = expand.grid(k=1:20),
      trControl = trControl,
      metric = "Accuracy",
      data = train[,-1])
saveRDS(tout, '../../Data/corp_knn.rds')
```

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

## Implementing KNN in R {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r}
tout
```

## KNN performance {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{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")
```

## Using KNN to fill in industry {data-background="../Backgrounds/case_digital.jpg" class="default"}

<div class="clearfix">
<div class="right">
```{r}
label$industry_pred <- predict(tout,
                              label)
label[,c("document",
         "industry_pred")] %>%
  head %>% html_df
```
</div>
<div class="left" style="font-size:24px">
1. <a href="https://www.sec.gov/Archives/edgar/data/817473/0000817473-14-000010-index.htm">American Capital</a>: Asset manager and private equity
    - SIC missing, but clearly finance <span class="hi1">✔</span>
2. <a href="https://www.sec.gov/Archives/edgar/data/52428/0000820027-14-000025-index.htm">Ameriprise Certificate Co</a>: Investment company
    - SIC missing, but clearly finance <span class="hi1">✔</span>
3. <a href="https://www.sec.gov/Archives/edgar/data/837465/0000837465-14-000002-index.htm">Callaway Golf</a>: Golf equipment
    - SIC 3949 <span class="hi1">✔</span>
4. <a href="https://www.sec.gov/Archives/edgar/data/837919/0000837919-14-000002-index.htm">Everest Fund L P</a>: Speculative trading of commodity futures
    - SIC 6221 <span class="hi1">✔</span>
5. <a href="https://www.sec.gov/Archives/edgar/data/9631/000089109214000570/0000891092-14-000570-index.htm">Bank of Nova Scotia</a>: Joint with Scotiabank Covered Bond Guarantor Limited Partnership
    - SIC 6022 <span class="hi1">✔</span>
6. <a href="https://www.sec.gov/Archives/edgar/data/1471824/0000891092-14-002078-index.htm">Teucrium Commodity Trust</a>: Commodity funds
    - SIC 6221 <span class="hi1">✔</span>
</div>
</div>


##  {data-background="../Backgrounds/case_digital.jpg" class="default"}

> "Any sufficiently advanced technology is indistinguishable from magic." -- Sir Arthur Charles Clarke

<img src="../Figures/arthur_clarke.jpg" height="400px">

## Bonus: t-SNE on KNN {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, echo=F, 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)
```

## Bonus: t-SNE on KNN (code) {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, eval=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)
}
```
```{r, eval=F}
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)
```



## Recap {data-background="../Backgrounds/ending.jpg" class="default"}

Today, we:

1. Processed a set of 6,000 annual reports from 2014 to examine their readability
2. Examined the content discussed in annual reports in 2014
3. Examined the natural groupings of content across firms
    - This doesn't necessarily match up well with SIC codes
    - There are some firms that don't quite fit with others in their industry (as we algorithmically identified)
4. Filled in missing industry data using KNN, and were correct in all 6 checked entries <span class="hi1">✔</span>

# End matter {data-background="../Backgrounds/ending.jpg" class="header"}

## For next week {data-background="../Backgrounds/ending.jpg" class="default"}

- For next week:
    - Datacamp
        - Do the assigned chapter on machine learning
    - Keep working on the group project
- Survey on the class session at <a href="https://rmc.link/420survey8">rmc.link/420survey8</a>

## Packages used for these slides {data-background="../Backgrounds/coding.jpg" class="default"}

- `package:caret`
- `package:cluster`
- `package:DT`
- `package:kableExtra`
- `package:knitr`
- `package:lattice`
- `package:quanteda` and `package:stopwords`
- `package:readtext`
- `package:revealjs`
- `package:Rtsne`
- `package:stm` and `package:stmBrowser`
- `package:tidyr`
- `package:tidyverse`
    - `package:dplyr`, `package:magrittr`, `package:readr`


## Custom code {data-background="../Backgrounds/coding.jpg" class="default"}

<small style="width:100%">

```{r,eval=F}
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)
  }
}
```
</small>



## Using more clusters {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, echo=F}
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()
```

## Using more clusters {data-background="../Backgrounds/coding.jpg" class="default"}

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

## Using more clusters {data-background="../Backgrounds/coding.jpg" class="default"}

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

## NOTES {class="never"}

- Simplify LDA slides -- too rambly, not simple enough
- Push the counts, dfm, and TF-IDF to be above LDA
    - Simple, then LDA.
- Add a simpler example, 2 dimensions, of kmeans
- Add a simpler example, 2 dimensions, of knn
    - Get rid of the kmeans slide on the knn section
- Add enron discussion

