---
title: "ACCT 420: Topic modeling and anomaly detection"
author: 
  - name: Dr. Richard M. Crowley
    email: rcrowley@smu.edu.sg
    affiliation: https://rmc.link/
title-slide-attributes:
  data-background-video: "../Backgrounds/title.mp4"
  data-background-video-loop: "true"
  data-background-size: stretch
format:
  revealjs:
    theme: default
    css: quarto_white.css
    navigation-mode: vertical
    controls: true
    transition: slide
    background-transition: fade
    transition-speed: default
    width: 1536
    height: 864
    logo: ../Logos/logo2.png
    code-link: true
    code-line-numbers: false
    chalkboard:
      theme: whiteboard
    auto-stretch: false
revealjs-plugins:
  - attribution
execute:
  echo: true
  warning: false
  message: false
editor: source
---

# Front Matter {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

<!--# Note that the downlit and xml2 pacakges are required to compile -->
<!--# Installed extensions include `attribution`, `qrcode`, and `material-icons` -->

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

<!--# List packages here that are not being correctly highlighted -->
<!--# The usual cause is calling the library from a echo=FALSE cell -->

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

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

:::: columns

::: {.column width="50%"}

![](../Course_Diagrams/Map-8.png){fig-alt="Roadmap for the course"}

:::

::: {.column width="50%"}

- **Theory:**
    - NLP
    - Anomaly detection
- **Application:**
    - Understand annual report readability
    - Examine the *content* of annual reports
    - Group firms on content
    - Fill in missing data
- **Methodology:**
    - ML/AI (LDA, k-means, KNN)
    - Dimensionality reduction: UMAP

:::

::::

## Group project tip \#1 {data-background="../Backgrounds/presentation.jpg"}

> For reading large files, `{readr}` is your friend

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

# OR

df <- read_csv("really_big_file.csv.gz")
```

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

## Group project tip \#2 {data-background="../Backgrounds/presentation.jpg"}

> For saving intermediary results, `saveRDS()` + `readRDS()`  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
    - As an added bonus, RDS files are compressed, taking less space on disk than csv files

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

# Sets of documents (corpus) {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

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

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

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

```{r, eval=FALSE, echo=FALSE, include=F}
temp <- summary(corp)
saveRDS(temp, "../../Data/Session_8-corp_summary.rds")
```

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

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

```{r, eval=FALSE, echo=FALSE, include=F}
temp <- textstat_readability(corp, "FOG")
saveRDS(temp, "../../Data/Session_8-corp_readability_FOG.rds")
```

```{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/Session_8-corp_readability_FOG.rds")
corp_FOG %>%
  head() %>%
  html_df()
```

> Recall that Micorsoft's annual report had a Fog index of 20.88

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

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

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

- 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: Microsoft is SIC 7372
    - 73: Business services
    - 737: Computer programming, data processing, and other computer related services
    - 7372: Prepackaged software

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

- Construct a data set of industries mapped to filings

<!--# Data is from WRDS SEC Analytics Suite - SEC Filings on WRDS -->

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

- Merge the industry data with the readability data

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

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

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

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

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

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

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



## `quanteda` bonus: References across text (Global warming) {data-background="../Backgrounds/default.jpg" .small}

```{r, eval=FALSE, echo=FALSE, include=F}
temp <- tokens(corp)
saveRDS(temp, '../../Data/Session_8-corp_tokens.rds')
temp2 <- kwic(temp, pattern = phrase("global warming"), window = 3)
saveRDS(temp2, "../../Data/Session_8-corp_kwic_gw.rds")
temp2 <- kwic(temp, pattern = phrase(c("COVID-19", "coronavirus")), window = 3)
saveRDS(temp2, "../../Data/Session_8-corp_kwic_covid.rds")
```


```{r, eval=F}
corp_tokens <- tokens(corp)  # This takes a couple hours to run

# kwic() is very fast to run though
kwic(corp_tokens, pattern = phrase("global warming"), window = 3) %>%
  as.tibble() %>%
  mutate(text=paste(pre,keyword,post)) %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  select(docname, text) %>%
  sample_n(100) %>%
  datatable(options = list(pageLength = 5), rownames=F)
```
```{r, echo=F, 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)
```

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

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

## `quanteda` bonus: References across text (COVID-19) {data-background="../Backgrounds/default.jpg" .small}

```{r, eval=F}
corp_tokens <- tokens(corp)  # This takes a couple hours to run

# kwic() is very fast to run though
kwic(corp_tokens, pattern = phrase(c("COVID-19", "coronavirus")), window = 3) %>%
  as.tibble() %>%
  mutate(text=paste(pre,keyword,post)) %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  select(docname, text) %>%
  sample_n(100) %>%
  datatable(options = list(pageLength = 5), rownames=F)
```


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

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

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

# Going beyond simple text measures {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

- 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"}

```{r, eval=F, include=F, echo=F}
words = 0
for (i in 1:length(corp_tokens)) {
  words = words + length(corp_tokens[[i]])
}
words

words / 10000 / 60 / 24
words / 700 / 60 / 24
words / 250 / 60 / 24
```

- This is a hard question to answer -- our sample has 317,759,360 words in it!
    - 22.1 days for the "world's fastest reader", per <a href="https://skeptoid.com/episodes/4229">this source</a>
    - 315.2 days for a standard speed reader (<a href="http://psycnet.apa.org/record/1998-11174-004">700wpm</a>)
    - **882.7 days** for an average reader (<a href="http://psycnet.apa.org/record/1998-11174-004">250wpm</a>)

::: {.callout-tip}
## Solutions

1. We could read a small sample of them.
    - Imprecise, risks missing out on some types of discussion
    - We need a second computer process to apply our findings to the rest of the documents
2. Have a computer read all of them!

:::


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

- 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"}

- 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"}

- In `{quanteda}`, use `dfm()`
- Useful additions:
    - We can pipe the output of `dfm()` to `dfm_remove()` to remove stopwords
        - You can use `remove=stopwords()` for a simple list
        - We can use SMART like last week: `remove=stopwords(source='smart')`
        - The `stopwords()` function is provided by the `{stopwords}` package, and actually supports over 50 languages, including Chinese, English, Hindi, and Malay
        - For other languages: `remove=stopwords("zh", source="stopwords-iso")`
        - With `remove=c(...)`, You can supply a list of stop words to remove
    - We can remove particularly frequent or infrequent terms with `dfm_trim()`
- We can preprocess our `tokens()` output as well
    - Pass it to `tokens_wordstem()` for stemming
        - Ex.: *cod*e, *cod*ing, and *cod*er would all become *cod*
    - `tokens()` has the options `remove_punct=T` and `remove_numbers=T` too

            

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

```{r, eval=FALSE, echo=FALSE, include=F}
corp_tokens2 <- tokens(corp_tokens, remove_punct=TRUE, remove_numbers=TRUE) %>%
  tokens_wordstem()

tdm <- dfm(corp_tokens2) %>%
  dfm_remove(stopwords(source='smart')) %>%
  dfm_trim(min_termfreq=10, termfreq_type = "count")

# adding industry to the tdm
docs <- docnames(corp)
docs <- data.frame(document=docs)
docs <- docs %>% left_join(df_SIC)
docvars(tdm, field="industry") <- docs$industry

temp <- topfeatures(tdm, n=5, groups=industry)
saveRDS(temp, '../../Data/Session_8-corp_dfm_feat.rds')

```

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

# With stopwords
tdm <- dfm(corp_tokens) %>%
       dfm_remove(stopwords(source='smart'))

# With stopwords and stemming -> Used in next slides
# 683M elements in the output
corp_tokens2 <- tokens(corp_tokens, remove_punct=TRUE, remove_numbers=TRUE) %>%
  tokens_wordstem()
tdm <- dfm(corp_tokens2) %>%
  dfm_remove(stopwords(source='smart'))
  dfm_trim(min_termfreq=10, termfreq_type = "count")
```

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



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

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

> This isn't very informative

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

- 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 `{quanteda}`
    - The options we'll specify are used to match a more standard output

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

$$
 \frac{f_{w,d}}{f_d} \cdot -\log_2\left(\frac{n_w}{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"}

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

saveRDS(temp, '../../Data/Session_8-corp_tfidf_feat.rds')
```

```{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/Session_8-corp_tfidf_feat.rds')
```

> These terms are often more meaningful

# Moving on to LDA {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

- **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

![](../Figures/dirichlet.png){height="100px" fig-align="center"}

> [More details from the creator](http://www.cs.columbia.edu/~blei/papers/Blei2012.pdf)

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

![](../Figures/LDA_blei.png){height="600px" fig-align="center"}

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

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 (MCMC method)

- It's a bit complicated mathematically, 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

::: {.callout-note}
## What type of Algorithm is LDA?

- Because of the distributional assumptions (which include priors), this is Bayesian
- Because of the way a Gibbs sampler approximates the distributions, this is machine learning

:::
   
## Implementations in R {data-background="../Backgrounds/reference.jpg"}

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



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

```{r, eval=FALSE, echo=FALSE, include=F}
out <- convert(tdm, to = 'stm')
saveRDS(out, '../../Data/Session_8-corp_out_stm.rds')
```

```{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/Session_8-corp_out_stm.rds')
#out <- readRDS('../../Data/Session_8-corp_out_lda.rds')
#out <- readRDS('../../Data/Session_8-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"}

- We will use the `stm()` function from the `{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 of LDA 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=FALSE, echo=FALSE, include=F}
library(stm)
topics <- stm(out$documents, out$vocab, K=10)
saveRDS(topics, '../../Data/Session_8-corp_stm_topics.rds')
```

```{r, eval=F}
library(stm)
topics <- stm(out$documents, out$vocab, K=10)
```
```{r, include=F}
library(stm)
topics <- readRDS('../../Data/Session_8-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"}

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

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

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

```{r, eval=FALSE, echo=FALSE, include=F}
toLDAvis(topics, out$documents, R=10, out.dir='LDAvis')

library(stmBrowser)
stmBrowser(topics, data=data.frame(text=names(out$documents),
                                   industry=out$meta$industry),
           c('industry'), text='text')
```

- Using LDAvis via `{STM}`'s `toLDAvis()` function
    - Need `{LDAvis}` and `{servr}` installed to run
    - Note: LDAvis scrambles the topic numbers (e.g., topic 1 is LDAvis' topic 9)

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

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

- Using `{stmBrowser}`'s `stmBrowser()` function
    - Install from github, not CRAN

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

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

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

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

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

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

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



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

- 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 2021 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"}

> 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-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

- 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"}

- 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

:::: {.columns}

::: {.column width="50%"}

- 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!

:::

::: {.column width="50%"}
![](../Figures/grouping_anom.png){height="250px" fig-align="center"}
:::

::::


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

$$
\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

:::: {.columns}

::: {.column width="50%"}

- Pros:
    - Very fast to run
    - Simple interpretation

:::

::: {.column width="50%"}

- Cons
    - Simple algorithm
    - Need to specify $k$, the number of clusters

:::

::::

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

- 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 `{tidyr}`


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

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

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

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

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
- Note: `kmeans()` is built into R -- no packages needed

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

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


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

- There is a relatively new method (2018), UMAP, that is significantly better
    - UMAP stands for **U**niform **M**anifold **A**pproximation and **P**rojection for Dimension Reduction
    - We will use it to reduce 68 dimensions down to 2
    - It is useful for plotting 2 dimensional representations of high dimensional data by maintaining *local* distance structures
        - It also maintains distances *globally*, mostly
    - It is computationally efficient
    - It is based on solid mathematical theory
        - Reimannian manifolds and geodesic distance

::: aside

There is also t-SNE (**t**-distributed **S**tochastic **N**eighbor **E**mbedding) from 2008, but it is inferior for 2 reasons: 1) it is more computationally costly than UMAP and 2) it is a bit misleading, as it only maintains distance locally *but not globally*. There is an even more outdated method (PCA), which struggles on higher dimensional data like our 10 topics.

:::
    
## Implementing UMAP {data-background="../Backgrounds/default.jpg"}

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

- We will use the `{uwot}` package to implement UMAP
- Our goal is to extract the locations that each document should be placed at in a 2D space
- The `umap()` function builds the model
- The `umap_train$embedding` object contains the needed coordinates
- Then we just add these back into our data

## Visualizing with UMAP: k-means {data-background="../Backgrounds/case_digital.jpg"}

:::: {.columns}

::: {.column width="50%"}

**Colored by SIC codes**

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

:::

::: {.column width="50%"}

**Colored by kmeans**

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

:::

::::

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

- 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"}

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

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

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

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

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

# Looking for anomalies {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

- 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(abs(mat - fitted(clusters))))  # Distance from center
wide_topics[,c(1,2,4,8,16)] %>% arrange(desc(dist)) %>% slice(1:5) %>% html_df()
```

> We are a blank check company incorporated on ________ as a Cayman Islands exempted company for the purpose of effecting a merger, share exchange, asset acquisition, share purchase, reorganization or similar business combination with one or more businesses or entities (a "Business Combination"). <cite>All 5 files...</cite>

- They are used for SPACs (e.g., Grab)

## Looking for anomalies (ignoring finance firms) {data-background="../Backgrounds/case_digital.jpg" .small}

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

- All: Yet more SPACs, just with the wrong industry in their filings...

- How many SPACs are there?

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

## Looking for anomalies (ignoring high compensation discussion) {data-background="../Backgrounds/case_digital.jpg" .small}

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

```

- 1: Sustainable homebuilder
- 2: Largest US homebuilder (4-7 are similar companies)
- 3: A bankrupt, regional lessor of 12 aircraft
- 8: Contracting services for automotive and energy firms; data center operation
- 9: A timeshare firm spun off from Hilton
- 10: A complex IPO-related entity with no actual operations

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

- 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>

## Brainstorming... {data-background="../Backgrounds/rain.jpg" class="default present-only"}

<iframe allowfullscreen frameborder="0" height="600px" mozallowfullscreen src="https://app.wooclap.com/events/KOJCMP/questions/5f85761b85d5387f093843e0" width="100%"></iframe>

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

> 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-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

- 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"}

- 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"}

- 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"}

- We'll use the `{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 %>% filter(!is.na(industry))
label <- wide_topics %>% filter(is.na(industry))
```

```{r, eval=F, echo=F, include=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[,c(2:11)])
saveRDS(tout, '../../Data/Session_8-corp_knn.rds')
```

```{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[,c(2:11)])
```

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

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

```{r}
tout
```

## KNN performance as we increase `k` {data-background="../Backgrounds/case_digital.jpg"}

```{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" .small}

:::: {.columns}

::: {.column width="50%"}

1. <a href="https://www.sec.gov/Archives/edgar/data/17313/000001731321000075/0000017313-21-000075-index.html">CAPITAL SOUTHWEST CORP</a>: "closed-end, non-diversified investment company"
    - 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">Rayonier Inc</a>: It is a timberland REIT, but it used to be a paper manufacturer
    - SIC is 6798 (finnace) for 1 entity, missing for another, but clearly finance <span class="hi2">x</span>
3. <a href="https://www.sec.gov/Archives/edgar/data/52428/000082002721000014/0000820027-21-000014-index.htm">AMERIPRISE CERTIFICATE COMPANY</a>: Financial certificate firm
    - SIC missing, but clearly finance <span class="hi1">✔</span>
4. <a href="https://www.sec.gov/Archives/edgar/data/837465/000083746521000003/0000837465-21-000003-index.htm">Callaway Golf</a>: Golf equipment
    - SIC 3949 (in manufacturing) <span class="hi1">✔</span>
5. Quest Management, Inc.: No operations, but used to do marketing for fitness equipment
    - No SIC, but it would fall under services <span class="hi1">✔</span>
6. <a href="https://www.sec.gov/Archives/edgar/data/1535778/000104746921000783/0001047469-21-000783-index.htm">MSC Income Fund</a>: "Closed end management investment company"
    - SIC missing, but clearly finance <span class="hi1">✔</span>

:::

::: {.column width="50%"}

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

:::

::::


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

Today, we:

1. Processed a set of 6,933 annual reports from 2021 to examine their readability
2. Examined the content discussed in annual reports in 2021
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 5 of 6 checked entries <span class="hi1">✔</span>







# End Matter {data-background-video="../Backgrounds/section.mp4"}

## Wrap up {data-background="../Backgrounds/ending.jpg"}

- Keep working on the project -- you have a lot of tools you can use already
    - And you will learn 1 more next week!
- Survey on the class session at this QR code:

{{< qrcode https://rmc.link/420survey8 qrsurvey width=200 height=200 >}}

## Packages used for these slides {data-background="../Backgrounds/code.jpg" visibility="uncounted"}

- `{caret}`
- `{cluster}`
- `{DT}`
- `{downlit}`
- `{kableExtra}`
- `{knitr}`
- `{quanteda}` and `{stopwords}`
- `{quarto}`
- `{readtext}`
- `{revealjs}`
- `{stm}` and `{stmBrowser}`
- `{tidyr}`
- `{tidyverse}`
- `{uwot}`


## NOTES {class="never" visibility="hidden"}

- Add a simpler example, 2 dimensions, of knn
- Add enron discussion
