---
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: 1536
      height: 864
      slideNumber: true
      controlsTutorial: true
      chalkboard:
        theme: whiteboard
        readOnly: false
---

## ACCT 420: Logistic Regression for Corporate Fraud {data-background-video="../Backgrounds/title.mp4" data-background-video-loop=true class="header no-toc-progress present-only" data-state="no-toc-progress"}

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

<h3 class="date"><br>Session 6</h3>
<h2 class="author"><br>Dr. Richard M. Crowley<br>rcrowley@smu.edu.sg<br><a href="http://rmc.link/">http://rmc.link/</a></h2>

## ACCT 420: Logistic Regression for Corporate Fraud {data-background="../Backgrounds/title.jpg" class="header no-toc-progress present-not" data-state="no-toc-progress"}

<h3 class="date"><br>Session 6</h3>
<h2 class="author"><br>Dr. Richard M. Crowley<br>rcrowley@smu.edu.sg<br><a href="http://rmc.link/">http://rmc.link/</a></h2>

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

```{r, include=F, eval=F}
# Run to prep data
df <- read.csv('../../Data/10-Ks_5_60-0-1100-100_LDA31_dechow_2004.csv')
df2 <- read.csv('../../Data/10-Ks_5_60-0-1100-100_LDA31_fin_2004.csv')
df3 <- left_join(df, df2[c(3:17,19,20,22:49)])  # col 49 is "Filing", which is a unique key
write.csv(df3, "../../Data/Session_6.csv", row.names=F)
df3_clean <- df3[,c(-125,-127)]
write.csv(df3_clean, "../../Data/Session_6-cleaned.csv", row.names=F)
```

```{r, include=F}
df <- read.csv("../../Data/Session_6.csv") %>%
  mutate(AAER = factor(AAER, levels=c(0,1)))
```

# 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-6.png" width="500px">
</div>
<div class="right objectives">
- <b>Theory:</b>
    - Economics
    - Psychology
- <b>Application:</b>
    - Predicting fraud contained in annual reports
- <b>Methodology:</b>
    - Logistic regression
    - LASSO
</div>
<div class="clearfix">

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

> - Explore on your own
> - No specific required class this week

# Corporate/Securities Fraud {data-background="../Backgrounds/section.jpg" class="header"}

## Traditional accounting fraud {data-background="../Backgrounds/fraud.jpg" class="default"}

1. A company is underperforming
2. Management cooks up some scheme to increase earnings
    - Worldcom (1999-2001)
        - Fake revenue entries
        - Capitalizing line costs (should be expensed)
    - Olympus (late 1980s-2011): Hide losses in a separate entity
        - "Tobashi scheme"
    - Wells Fargo (2011-2018?)
        - Fake/duplicate customers and transactions
3. Create accounting statements using the fake information

## Reversing it {data-background="../Backgrounds/cookie_jar.jpg" class="default"}

1. A company is overperforming
2. Management cooks up a scheme to "save up" excess performance for a rainy day
    - <a href="https://www.economist.com/newsbook/2010/07/23/taking-away-dells-cookie-jar">Dell (2002-2007)</a>
        - Cookie jar reserve, from secret payments by Intel, made up to **76%** of quarterly income
    - <a href="https://www.sec.gov/news/press/2004-105.htm">Brystol-Myers Squibb (2000-2001)</a>
3. Recognize revenue/earnings when needed in the future to hit earnings targets


## Other accounting fraud types {data-background="../Backgrounds/fraud.jpg" class="default"}

- *Options backdating*: <a href="https://www.sec.gov/news/press/2007/2007-70.htm"> Apple (2001)</a>
- Using an auditor that *isn't registered*: <a href="https://dart.deloitte.com/USDART/resource/b44c3afb-3f7f-11e6-95db-51a9f8be3f47">Commerce Group Corp (2003)</a>
- Releasing financial statements that were *not reviewed by an auditor*: <a href="https://www.sec.gov/litigation/admin/2018/34-84258.pdf">Cardiff International (2017)</a>
- *Related party transactions* (transferring funds to family members): <a href="https://www.sec.gov/litigation/litreleases/2012/lr22552.htm">China North East Petroleum Holdings Limited</a>
- *Insufficient internal controls*: <a href="https://www.sec.gov/litigation/admin/2018/34-83858.pdf">Citigroup (2008-2014)</a> via Banamex and <a href="http://eresources.nlb.gov.sg/infopedia/articles/SIP_422_2005-01-25.html">Asia Pacific Breweries</a>
- *Round-tripping*: Transactions to inflate revenue that have no substance: <a href="https://www.sec.gov/news/press/2004-2.htm">Suprema Specialties (1998-2001)</a>
- *Bribery*: <a href="https://www.channelnewsasia.com/news/business/keppel-o-m-bribery-case-what-you-need-to-know-9836154">Keppel O&M (2001-2014)</a>, $55M USD in bribes to Brazilian officials for contracts
- *Fake the whole company*: <a href="https://www.nytimes.com/1990/02/25/books/nothing-but-zzzz-best.html">ZZZZ Best (1982-1987)</a>
    - Getting funding from insurance fraud, theft, credit card fraud, and fake contracts; faking a real project to get a clean audit to take the company public
- *Ponzi scheme*: <a href="https://www.nytimes.com/2009/01/25/business/25bernie.html">Bernard Madoff</a>
- *Material omissions and misstatements*: <a href="https://www.sec.gov/litigation/litreleases/2013/lr22801.htm">Imaging Diagnostic Systems (2013)</a>
- *Failed to file annual and quarterly reports*: <a href="https://www.sec.gov/litigation/admin/2010/34-61344a.pdf">Applied Wellness Corporation (2008)</a>
- *Aiding another company's fraud* (Take Two, by parking 2 video games): <a href="https://www.sec.gov/litigation/admin/2008/34-57303.pdf">Capitol Distributing LLC</a>
- *Misleading statements on Twitter*: <a href="https://www.sec.gov/news/press-release/2018-219">Tesla (2018)</a>

## Some of the more interesting cases {data-background="../Backgrounds/fraud.jpg" class="default"}

- <a href="https://www.sec.gov/litigation/admin/3437730.txt">AMD (1992-1993)</a>
    - Claimed it was developing processor microcode independently, when it actually provided Intel's microcode to it's engineers
- <a href="https://www.sec.gov/litigation/litreleases/lr17024.htm">Am-Pac International (1997)</a>
    - Sham sale-leaseback of a bar to a corporate officer
- <a href="https://www.sec.gov/litigation/admin/2007/33-8815.pdf">CVS (2000)</a>
    - Not using mark-to-market accounting to fair value stuffed animal inventories
- <a href="https://www.sec.gov/litigation/litreleases/lr16732.htm">Countryland Wellness Resorts, Inc. (1997-2000)</a>
    - Gold reserves were actually... dirt.
- <a href="https://www.straitstimes.com/singapore/courts-crime/keppel-club-duo-convicted-for-37m-membership-scam">Keppel Club (2014)</a>
    - Employees created 1,280 fake memberships, sold them, and retained all profits (\$37.5M)
    
## What will we look at today? {data-background="../Backgrounds/default.jpg" class="default"}

> Misstatements: Errors that affect firms' accounting statements or disclosures which were done seemingly *intentionally* by management or other employees at the firm.

<img src="../Figures/johnny-automatic-bag-of-money.png" height="150px">

## How do misstatements come to light? {data-background="../Backgrounds/default.jpg" class="default"}

1. The company/management admits to it publicly
2. A government entity forces the company to disclose
    - In more egregious cases, government agencies may disclose the fraud publicly as well
3. Investors sue the firm, forcing disclosure

## Where are these disclosed? (US) {data-background="../Backgrounds/default.jpg" class="default"}

1. <a href="https://www.sec.gov/divisions/enforce/friactions.shtml">US SEC AAERs</a>: Accounting and Auditing Enforcement Releases
    - Highlight larger/more important cases, written by the SEC
    - Example: The *Summary* section of <a href="https://www.sec.gov/litigation/admin/2018/34-84017.pdf">this AAER against Sanofi</a>
2. 10-K/A filings ("10-K" $\Rightarrow$ annual report, "/A" $\Rightarrow$ amendment)
    - Note: not all 10-K/A filings are caused by fraud!
        - Benign corrections or adjustments can also be filed as a 10-K/A
        - Note: <a href="https://www.auditanalytics.com/blog/reasons-for-an-amended-10-k-2017/">Audit Analytics' write-up on this for 2017</a>
3. By the US government through a 13(b) action
4. In a note inside a 10-K filing
    - These are sometimes referred to as "little r" restatements
5. In a press release, which is later filed with the US SEC as an 8-K
    - 8-Ks are filed for many other reasons too though

## Where are we at? {data-background="../Backgrounds/default.jpg" class="default"}

> Fraud happens in many ways, for many reasons

- All of them are important to capture
- All of them affect accounting numbers differently
- None of the individual methods are frequent...

> It is disclosed in many places.  All have subtly different meanings and implications

- We need to be careful here (or check multiple sources)

> This is a hard problem!

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

- Today we will examine these AAERs
    - Using a proprietary data set of >1,000 such releases
- To get a sense of the data we're working with, read the *Summary* section (starting on page 2) of this AAER against Sanofi
    - <a href="https://www.sec.gov/litigation/admin/2018/34-84017.pdf">rmc.link/420class6</a>

> Why did the SEC release this AAER regarding Sanofi?

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

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

# Predicting Fraud {data-background="../Backgrounds/section.jpg" class="header"}

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

> How can we *detect* if a firm *is* involved in a major instance of missreporting?

- This is a pure forensic analytics question
- "Major instance of misreporting" will be implemented using AAERs

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

- In these slides, I'll walk through the primary detection methods since the 1990s, up to currently used methods
- 1990s: Financials and financial ratios
    - Follow up in 2011
- Late 2000s/early 2010s: Characteristics of firm's disclosures
- mid 2010s: More holistic text-based measures of disclosures
    - This will tie to next lesson where we will explore how to work with text

> All of these are discussed in a <a href="https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2803733">Brown, Crowley and Elliott (2020 JAR)</a> -- I will refer to the paper as **BCE** for short

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

- I have provided some preprocessed data, sanitized of AAER data (which is partially public, partially proprietary)
- It contains 401 variables
     - From Compustat, CRSP, and the SEC (which I personally collected)
     - Many precalculated measures including:
         - Firm characteristics, such as auditor type (`bigNaudit`, `midNaudit`)
         - Financial measures, such as total accruals (`rsst_acc`)
         - Financial ratios, such as ROA (`ni_at`)
         - Annual report characteristics, such as the mean sentence length (`sentlen_u`)
         - Machine learning based content analysis (everything with `Topic_` prepended)

> Pulled from BCE's working files

## Training and Testing {data-background="../Backgrounds/question.jpg" class="default"}

- Already has testing and training set up in variable `Test`
    - Training is annual reports released in 1999 through 2003
    - Testing is annual reports released in 2004

> What potential issues are there with our usual training and testing strategy?

<div class="fragment exclude">

- There is a significant lag between when a fraud is caught and when a fraud actually happened!
    - To mirror the available information in 2004, we should *censor* `AAER` for training such that it only captures AAERs known by 2003

```{r, echo=F}
ex <- data.frame(year=c(1999,2001,2003), year_found=c(2001,2003,2006), aaer=c(1,1,1), aaer_as_of_2004=c(1,1,0))
html_df(ex)
```

</div>

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

- Censoring training data helps to emulate historical situations
    - Build an algorithm using only the data that was available at the time a decision would need to have been made
- **Do not censor the testing data**
    - Testing emulates where we want to make an optimal choice in real life
        - We want to find frauds regardless of how well hidden they are!

<img src="../Figures/Censoring.png" width="800px">

## Event frequency {data-background="../Backgrounds/question.jpg" class="default"}

- Very low event frequencies can make things tricky

```{r}
df %>%
  group_by(year) %>%
  mutate(total_AAERS = sum(AAER==1), total_observations=n()) %>%
  slice(1) %>%
  ungroup() %>%
  select(year, total_AAERS, total_observations) %>%
  html_df
```

> 246 AAERs in the training data, 401 total variables...

## Dealing with infrequent events {data-background="../Backgrounds/reference.jpg" class="default"}

- A few ways to handle this
    1. Very careful model selection (keep it sufficiently simple)
    2. Sophisticated degenerate variable identification criterion + simulation to implement complex models that are just barely simple enough
        - The main method in BCE
    3. Automated methodologies for pairing down models
        - We'll discuss using LASSO for this at the end of class
            - Also implemented in BCE

# 1990s approach {data-background="../Backgrounds/section.jpg" class="header"}

## The 1990s model {data-background="../Backgrounds/case.jpg" class="default"}

- Many financial measures and ratios can help to predict fraud

<div class="clearfix">
<div class="left">
- EBIT
- Earnings / revenue
- ROA
- Log of liabilities
- liabilities / equity
- liabilities / assets
- quick ratio
- Working capital / assets
- Inventory / revenue
- inventory / assets
- earnings / PP&E
- A/R / revenue
</div>
<div class="right">
- Change in revenue
- Change in A/R + 1
- $>10\%$ change in A/R
- Change in gross profit + 1
- $>10\%$ change in gross profit
- Gross profit / assets
- Revenue minus gross profit
- Cash / assets
- Log of assets
- PP&E / assets
- Working capital
</div>
</div>

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

- Purely economic
- Misreporting firms' financials should be different than expected
    - Perhaps more income
    - Odd capital structure
    - Odd balance of receivables

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

```{r, warning=F, }
fit_1990s <- glm(AAER ~ ebit + ni_revt + ni_at + log_lt + ltl_at + lt_seq +
                   lt_at + act_lct + aq_lct + wcap_at + invt_revt + invt_at +
                   ni_ppent + rect_revt + revt_at + d_revt + b_rect + b_rect +
                   r_gp + b_gp + gp_at + revt_m_gp + ch_at + log_at +
                   ppent_at + wcap,
                 data=df[df$Test==0,],
                 family=binomial)
summary(fit_1990s)
```

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

```{r, echo=F, warning=F, message=F, fig.height=5}
library(yardstick)

df$pred_1990s <- predict(fit_1990s, df, type="response")

auc_in <- df %>% filter(Test==0) %>% roc_auc(truth=AAER, estimate=pred_1990s, event_level='second')
auc_out <- df %>% filter(Test==1) %>% roc_auc(truth=AAER, estimate=pred_1990s, event_level='second')
curve_in_1990s <- df %>% filter(Test==0) %>% roc_curve(truth=AAER, estimate=pred_1990s, event_level='second')
curve_out_1990s <- df %>% filter(Test==1) %>% roc_curve(truth=AAER, estimate=pred_1990s, event_level='second')

ggplot() +
  geom_line(data=curve_in_1990s, aes(y=sensitivity, x=1-specificity, color="1990s, In Sample")) + 
  geom_line(data=curve_out_1990s, aes(y=sensitivity, x=1-specificity, color="1990s, Out of Sample")) +
  geom_abline(slope=1)

aucs_1990s <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_1990s) <- c("In sample AUC", "Out of sample AUC")
aucs_1990s
```

# The 2011 follow up {data-background="../Backgrounds/section.jpg" class="header"}

## The 2011 model {data-background="../Backgrounds/case.jpg" class="default"}

<div class="clearfix">
<div class="left">
- Log of assets
- Total accruals
- % change in A/R
- % change in inventory
- % soft assets
- % change in sales from cash
- % change in ROA
- Indicator for stock/bond issuance
- Indicator for operating leases
- BV equity / MV equity
</div>
<div class="right">
- Lag of stock return minus value weighted market return
- **Below are BCE's additions**
- Indicator for mergers
- Indicator for Big N auditor
- Indicator for medium size auditor
- Total financing raised
- Net amount of new capital raised 
- Indicator for restructuring
</div>
</div>

> Based on <a href="https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1911-3846.2010.01041.x">Dechow, Ge, Larson and Sloan (2011)</a>

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

```{r, warning=F}
fit_2011 <- glm(AAER ~ logtotasset + rsst_acc + chg_recv + chg_inv +
                  soft_assets + pct_chg_cashsales + chg_roa + issuance +
                  oplease_dum + book_mkt + lag_sdvol + merger + bigNaudit +
                  midNaudit + cffin + exfin + restruct,
                 data=df[df$Test==0,],
                 family=binomial)
summary(fit_2011)
```

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

```{r, echo=F, warning=F, message=F, fig.height=5}
df$pred_2011 <- predict(fit_2011, df, type="response")

auc_in <- df %>% filter(Test==0) %>% roc_auc(truth=AAER, estimate=pred_2011, event_level='second')
auc_out <- df %>% filter(Test==1) %>% roc_auc(truth=AAER, estimate=pred_2011, event_level='second')
curve_in_2011 <- df %>% filter(Test==0) %>% roc_curve(truth=AAER, estimate=pred_2011, event_level='second')
curve_out_2011 <- df %>% filter(Test==1) %>% roc_curve(truth=AAER, estimate=pred_2011, event_level='second')

ggplot() +
  geom_line(data=curve_in_2011, aes(y=sensitivity, x=1-specificity, color="2011, In Sample")) + 
  geom_line(data=curve_out_2011, aes(y=sensitivity, x=1-specificity, color="2011, Out of Sample")) +
  geom_abline(slope=1)

aucs_2011 <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_2011) <- c("In sample AUC", "Out of sample AUC")
aucs_2011
```

```{r, include=F, echo=F, warning=F, message=F}
df_save <- data.frame(Test=df$Test, AAER=df$AAER, pred_F=df$pred_2011)
```

# Late 2000s/early 2010s approach {data-background="../Backgrounds/section.jpg" class="header"}

## The late 2000s/early 2010s model {data-background="../Backgrounds/case.jpg" class="default"}

<div class="clearfix">
<div class="left">
- Log of # of bullet points + 1
- \# of characters in file header
- \# of excess newlines
- Amount of html tags
- Length of cleaned file, characters
- Mean sentence length, words
- S.D. of word length
- S.D. of paragraph length (sentences)
</div>
<div class="right">
- Word choice variation
- Readability
    - Coleman Liau Index
    - Fog Index
- % active voice sentences
- % passive voice sentences
- \# of all cap words
- \# of !
- \# of ?
</div>
</div>

> From a variety of papers

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

- Generally pulled from the communications literature
    - Sometimes ad hoc
- The main idea:
    - Companies that are misreporting probably write their annual report differently

## The late 2000s/early 2010s model {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, warning=F}
fit_2000s <- glm(AAER ~ bullets + headerlen + newlines + alltags +
                   processedsize + sentlen_u + wordlen_s + paralen_s +
                   repetitious_p + sentlen_s + typetoken + clindex + fog +
                   active_p + passive_p + lm_negative_p + lm_positive_p +
                   allcaps + exclamationpoints + questionmarks,
                 data=df[df$Test==0,],
                 family=binomial)
summary(fit_2000s)
```

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

```{r, echo=F, warning=F, message=F, fig.height=5}
df$pred_2000s <- predict(fit_2000s, df, type="response")

auc_in <- df %>% filter(Test==0) %>% roc_auc(truth=AAER, estimate=pred_2000s, event_level='second')
auc_out <- df %>% filter(Test==1) %>% roc_auc(truth=AAER, estimate=pred_2000s, event_level='second')
curve_in_2000s <- df %>% filter(Test==0) %>% roc_curve(truth=AAER, estimate=pred_2000s, event_level='second')
curve_out_2000s <- df %>% filter(Test==1) %>% roc_curve(truth=AAER, estimate=pred_2000s, event_level='second')

ggplot() +
  geom_line(data=curve_in_2000s, aes(y=sensitivity, x=1-specificity, color="2000s, In Sample")) + 
  geom_line(data=curve_out_2000s, aes(y=sensitivity, x=1-specificity, color="2000s, Out of Sample")) +
  geom_abline(slope=1)

aucs_2000s <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_2000s) <- c("In sample AUC", "Out of sample AUC")
aucs_2000s
```

```{r, include=F, echo=F, warning=F, message=F}
df_save$pred_S <- df$pred_2000s
```

## Combining the 2000s and 2011 models  {data-background="../Backgrounds/question.jpg" class="default"}

> Why is it appropriate to combine the 2011 model with the 2000s model?

- 2011 model: Parsimonious financial model
- 2000s model: Textual characteristics

<div class="fragment exclude">

> Little theoretical overlap

> Limited multicollinearity across measures

</div>

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

```{r, warning=F}
fit_2000f <- glm(AAER ~ logtotasset + rsst_acc + chg_recv + chg_inv +
                   soft_assets + pct_chg_cashsales + chg_roa + issuance +
                   oplease_dum + book_mkt + lag_sdvol + merger + bigNaudit +
                   midNaudit + cffin + exfin + restruct + bullets + headerlen +
                   newlines + alltags + processedsize + sentlen_u + wordlen_s +
                   paralen_s + repetitious_p + sentlen_s + typetoken +
                   clindex + fog + active_p + passive_p + lm_negative_p +
                   lm_positive_p + allcaps + exclamationpoints + questionmarks,
                 data=df[df$Test==0,],
                 family=binomial)
summary(fit_2000f)
```

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

```{r, echo=F, warning=F, message=F, fig.height=5}
df$pred_2000f <- predict(fit_2000f, df, type="response")

auc_in <- df %>% filter(Test==0) %>% roc_auc(truth=AAER, estimate=pred_2000f, event_level='second')
auc_out <- df %>% filter(Test==1) %>% roc_auc(truth=AAER, estimate=pred_2000f, event_level='second')
curve_in_2000f <- df %>% filter(Test==0) %>% roc_curve(truth=AAER, estimate=pred_2000f, event_level='second')
curve_out_2000f <- df %>% filter(Test==1) %>% roc_curve(truth=AAER, estimate=pred_2000f, event_level='second')

ggplot() +
  geom_line(data=curve_in_2000f, aes(y=sensitivity, x=1-specificity, color="2000s + 2011, In Sample")) + 
  geom_line(data=curve_out_2000f, aes(y=sensitivity, x=1-specificity, color="2000s + 2011, Out of Sample")) +
  geom_abline(slope=1)

aucs_2000f <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_2000f) <- c("In sample AUC", "Out of sample AUC")
aucs_2000f
```

```{r, include=F, echo=F, warning=F, message=F}
df_save$pred_FS <- df$pred_2000f
```

# The BCE model {data-background="../Backgrounds/section.jpg" class="header"}

## The BCE approach {data-background="../Figures/LDA2012.png" class="default"}

- Retain the variables from the other regressions
- Add in a machine-learning based measure quantifying how much documents talked about different topics common across all filings
    - Learned on just the 1999-2003 filings

## What the topics look like {data-background="../Backgrounds/default.jpg" class="default"}

<img src="../Figures/LDA2004.png" width="800px">

## Theory behind the BCE model {data-background="../Backgrounds/question.jpg" class="default"}

> Why use document content?

<div class="fragment exclude">

- From communications and psychology:
    - When people are trying to deceive others, what they say is carefully picked
        - Topics chosen are intentional
- Putting this in a business context:
    - If you are manipulating inventory, you don't talk about it

</div>

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

```{r, warning=F}
BCE_eq = as.formula(paste("AAER ~ logtotasset + rsst_acc + chg_recv + chg_inv +
  soft_assets + pct_chg_cashsales + chg_roa + issuance +
  oplease_dum + book_mkt + lag_sdvol + merger + bigNaudit +
  midNaudit + cffin + exfin + restruct + bullets + headerlen +
  newlines + alltags + processedsize + sentlen_u + wordlen_s +
  paralen_s + repetitious_p + sentlen_s + typetoken +
  clindex + fog + active_p + passive_p + lm_negative_p +
  lm_positive_p + allcaps + exclamationpoints + questionmarks + ",
  paste(paste0("Topic_",1:30,"_n_oI"), collapse=" + "), collapse=""))
fit_BCE <- glm(BCE_eq,
               data=df[df$Test==0,],
               family=binomial)
summary(fit_BCE)
```

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

```{r, echo=F, warning=F, message=F, fig.height=5}
df$pred_BCE <- predict(fit_BCE, df, type="response")

auc_in <- df %>% filter(Test==0) %>% roc_auc(truth=AAER, estimate=pred_BCE, event_level='second')
auc_out <- df %>% filter(Test==1) %>% roc_auc(truth=AAER, estimate=pred_BCE, event_level='second')
curve_in_BCE <- df %>% filter(Test==0) %>% roc_curve(truth=AAER, estimate=pred_BCE, event_level='second')
curve_out_BCE <- df %>% filter(Test==1) %>% roc_curve(truth=AAER, estimate=pred_BCE, event_level='second')

ggplot() +
  geom_line(data=curve_in_BCE, aes(y=sensitivity, x=1-specificity, color="BCE, In Sample")) + 
  geom_line(data=curve_out_BCE, aes(y=sensitivity, x=1-specificity, color="BCE, Out of Sample")) +
  geom_abline(slope=1)

aucs_BCE <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_BCE) <- c("In sample AUC", "Out of sample AUC")
aucs_BCE
```

```{r, include=F, echo=F, warning=F, message=F}
df_save$pred_BCE <- df$pred_aucs_BCE
```

## Comparison across all models {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, echo=F, warning=F, message=F, fig.height=5}

ggplot() +
  geom_line(data=curve_out_BCE, aes(y=sensitivity, x=1-specificity, color="BCE")) +
  geom_line(data=curve_out_2000f, aes(y=sensitivity, x=1-specificity, color="2000s + 2011")) +
  geom_line(data=curve_out_2000s, aes(y=sensitivity, x=1-specificity, color="2000s")) +
  geom_line(data=curve_out_2011, aes(y=sensitivity, x=1-specificity, color="2011")) +
  geom_line(data=curve_out_1990s, aes(y=sensitivity, x=1-specificity, color="1990s")) +
  geom_abline(slope=1) + 
  ggtitle("Out of Sample ROC Curves")

oos_aucs <- c(aucs_1990s[2], aucs_2011[2], aucs_2000s[2], aucs_2000f[2], aucs_BCE[2])
names(oos_aucs) <- c("1990s", "2011", "2000s", "2000s + 2011", "BCE")
oos_aucs
```


# Simplifying models with LASSO {data-background="../Backgrounds/section.jpg" class="header"}

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

- Least Absolute Shrinkage and Selection Operator
    - Least absolute: uses an error term like $|\varepsilon|$
    - Shrinkage: it will make coefficients smaller
        - Less sensitive → less overfitting issues
    - Selection: it will completely remove some variables
        - Less variables → less overfitting issues
- Sometimes called $L^1$ regularization
    - $L^1$ means 1 dimensional distance, i.e., $|\varepsilon|$

> Great if you have way too many inputs in your model

- This is how we can, in theory, put more variables in our model than data points

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

$$
\min_{\beta\in \mathbb{R}} \left\{\frac{1}{N}\left|\varepsilon\right|_2^2 + \lambda \left|\beta\right|_1\right\}
$$
<div class="clearfix">
<div class="left">
- Add an additional penalty term that is increasing in the absolute value of each $\beta$
    - Incentivizes lower $\beta$s, *shrinking* them
- The selection is part is explainable geometrically
</div>
<div class="right">
<img src="../Figures/LASSO.png" width="400px">
</div>
</div>

## Why use it? {data-background="../Backgrounds/reference.jpg" class="default"}

1. We have a preference for simpler models
2. Some problems are naturally very complex
    - Many linkages between different theoretical constructs
3. We don't have a good judgment on what theories are better than others for the problem

> LASSO lets us implement all of our ideas, and then it econometrically kicks out the ineffective ideas (model selection)

## Package for LASSO {data-background="../Backgrounds/default.jpg" class="default"}

- `package:glmnet`
  1. For all regression commands, they expect a `y` vector and an `x` matrix instead of our usual `y ~ x` formula
      - R has a helper function to convert a formula to a matrix: `model.matrix()`
          - Supply it the right hand side of the equation, starting with `~`, and your data
          - It outputs the matrix `x`
      - Alternatively, use `as.matrix()` on a data frame of your input variables
  2. It's family argument should be specified in quotes, i.e., `"binomial"` instead of `binomial`

## What else can the package do? {data-background="../Backgrounds/default.jpg" class="default"}

<div class="clearfix">
<div class="left">
Ridge regression

- Similar to LASSO, but with an $L^2$ penalty (Euclidean norm)

<img src="../Figures/ridge.png" width="400px">
</div>
<div class="right">
Elastic net regression

- Hybrid of LASSO and Ridge
- Below image by <a href="https://jaredlander.com/content/2015/11/LassoForEveryone.html">Jared Lander</a>

<img src="../Figures/elastic_net.png" width="400px">

</div>
</div>

## How to run a LASSO {data-background="../Backgrounds/case_digital.jpg" class="default"}

- To run a simple LASSO model, use `glmnet()`
- Let's LASSO the BCE model

```{r, message=F}
library(glmnet)
x <- model.matrix(BCE_eq, data=df[df$Test==0,])[,-1]  # [,-1] to remove intercept
y <- model.frame(BCE_eq, data=df[df$Test==0,])[,"AAER"]
fit_LASSO <- glmnet(x=x, y=y,
                    family = "binomial",
                    alpha = 1  # Specifies LASSO.  alpha = 0 is ridge
                    )
```

- Note: the model selection can be more elegantly done using the `package:useful` package, <a href="https://www.jaredlander.com/2018/02/using-coefplot-with-glmnet/">see here for an example</a>

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

```{r, fig.height=4}
plot(fit_LASSO)
```

## What's under the hood? {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r}
print(fit_LASSO)
```

## **One** of the 100 models {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, warning=F, fig.height=5}
#coef(fit_LASSO, s=0.002031)
coefplot(fit_LASSO, lambda=0.002031, sort='magnitude')
```

## How does this perform? {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r}
# na.pass has model.matrix retain NA values (so the # of rows is constant)
xp <- model.matrix(BCE_eq, data=df, na.action='na.pass')[,-1]
# s= specifies the version of the model to use
df$pred_L1 <- c(predict(fit_LASSO, xp, type="response", s = 0.002031))
```

```{r, echo=F, fig.height=4}
auc_in <- df %>% filter(Test==0) %>% roc_auc(truth=AAER, estimate=pred_L1, event_level='second')
auc_out <- df %>% filter(Test==1) %>% roc_auc(truth=AAER, estimate=pred_L1, event_level='second')
curve_in_L1 <- df %>% filter(Test==0) %>% roc_curve(truth=AAER, estimate=pred_L1, event_level='second')
curve_out_L1 <- df %>% filter(Test==1) %>% roc_curve(truth=AAER, estimate=pred_L1, event_level='second')

ggplot() +
  geom_line(data=curve_in_L1, aes(y=sensitivity, x=1-specificity, color="BCE, LASSO, In Sample")) + 
  geom_line(data=curve_out_L1, aes(y=sensitivity, x=1-specificity, color="BCE, LASSO, Out of Sample")) +
  geom_abline(slope=1)

aucs_L1 <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_L1) <- c("In sample AUC", "Out of sample AUC")
aucs_L1
```

## Automating model selection {data-background="../Backgrounds/reference.jpg" class="default"}

- LASSO seems nice, but picking between the 100 models is tough!
- It also contains a method of *$k$-fold cross validation* (default, $k=10$)
    1. Randomly splits the data into $k$ groups
    2. Runs the algorithm on 90% of the data ($k-1$ groups)
    3. Determines the best model
    4. Repeat steps 2 and 3 $k-1$ more times
    5. Uses the best overall model across all $k$ hold out samples
- It gives 2 model options:
    - `"lambda.min"`: The best performing model
    - `"lambda.1se"`: The simplest model within 1 standard error of `"lambda.min"`
        - This is the better choice if you are concerned about overfitting

## Running a cross validated model {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r}
# Cross validation
set.seed(697435)  #for reproducibility
cvfit = cv.glmnet(x=x, y=y,family = "binomial", alpha = 1, type.measure="auc")
```

<div class="clearfix">
<div class="left">
```{r}
plot(cvfit)
```
</div>
<div class="right">
```{r}
cvfit$lambda.min
cvfit$lambda.1se
```

> These are the dashed lines on the plot
</div>
</div>

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

<div class="clearfix">
<div class="left">
`lambda.min`
```{r, warning=F, fig.height=5.2, fig.width=4, echo=F}
#coef(cvfit, s = "lambda.min")
coefplot(cvfit, lambda='lambda.min', sort='magnitude') + theme(axis.text.y = element_text(size=15))
```
</div>
<div class="right">
`lambda.1se`
```{r, warning=F, fig.height=5.2, fig.width=4, echo=F}
#coef(cvfit, s = "lambda.1se")
coefplot(cvfit, lambda='lambda.1se', sort='magnitude') + theme(axis.text.y = element_text(size=15))
```
</div>
</div>

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

```{r}
# s= specifies the version of the model to use
df$pred_L1.min <- c(predict(cvfit, xp, type="response", s = "lambda.min"))
df$pred_L1.1se <- c(predict(cvfit, xp, type="response", s = "lambda.1se"))
```

```{r, include=F, echo=F, warning=F, message=F}
df_save$pred_lmin <- df$pred_L1.min
df_save$pred_l1se <- df$pred_L1.1se
```

```{r, echo=F, fig.height=3.5}

auc_in <- df %>% filter(Test==0) %>% roc_auc(truth=AAER, estimate=pred_L1.min, event_level='second')
auc_out <- df %>% filter(Test==1) %>% roc_auc(truth=AAER, estimate=pred_L1.min, event_level='second')
curve_in_CV.min <- df %>% filter(Test==0) %>% roc_curve(truth=AAER, estimate=pred_L1.min, event_level='second')
curve_out_CV.min <- df %>% filter(Test==1) %>% roc_curve(truth=AAER, estimate=pred_L1.min, event_level='second')

aucs_CV.min <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_CV.min) <- c("In sample, lambda.min", "Out of sample, lambda.min")

auc_in <- df %>% filter(Test==0) %>% roc_auc(truth=AAER, estimate=pred_L1.1se, event_level='second')
auc_out <- df %>% filter(Test==1) %>% roc_auc(truth=AAER, estimate=pred_L1.1se, event_level='second')
curve_in_CV.1se <- df %>% filter(Test==0) %>% roc_curve(truth=AAER, estimate=pred_L1.1se, event_level='second')
curve_out_CV.1se <- df %>% filter(Test==1) %>% roc_curve(truth=AAER, estimate=pred_L1.1se, event_level='second')

aucs_CV.1se <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_CV.1se) <- c("In sample, lambda.1se", "Out of sample, lambda.1se")

ggplot() +
  geom_line(data=curve_in_CV.min, aes(y=sensitivity, x=1-specificity, color="In Sample, lambda.min")) + 
  geom_line(data=curve_out_CV.min, aes(y=sensitivity, x=1-specificity, color="Out of Sample, lambda.min")) +
  geom_line(data=curve_in_CV.1se, aes(y=sensitivity, x=1-specificity, color="In Sample, lambda.1se")) + 
  geom_line(data=curve_out_CV.1se, aes(y=sensitivity, x=1-specificity, color="Out of Sample, lambda.1se")) +
  geom_abline(slope=1)

c(aucs_CV.min, aucs_CV.1se)
```

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

1. No p-values on coefficients
    - Simple solution -- run the resulting model with `glm()`
    - Solution only if using `family="gaussian"`:
        - Run the lasso use the `package:lars` package
            - `m <- lars(x=x, y=y, type="lasso")`
        - Then test coefficients using the `package:covTest` package
            - `covTest(m, x, y)`
2. Generally worse in sample performance
3. Sometimes worse out of sample performance (short run)
    - BUT: predictions will be more stable

# Wrap up {data-background="../Backgrounds/section.jpg" class="header"}

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

> What other data could we use to predict corporate fraud?

- What is the reason that this event or data would be useful for prediction?
    - I.e., how does it fit into your mental model?
- What if we were...
    - Auditors?
    - Internal auditors?
    - Regulators?
    - Investors?

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

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

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

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

- Next week:
    - Third assignment
        - On binary prediction
        - Finish in three weeks
        - Can be done in pairs
        - Submit on eLearn
    - Datacamp
        - Practice a bit more to keep up to date
            - Using R more will make it more natural
- Survey on the class session at <a href="https://rmc.link/420survey6">rmc.link/420survey6</a>

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

> Predicting class action lawsuits

- Another question that has both forecasting and forensic flair to it
    - Forensic: Often these companies were doing something wrong for a while in the past
    - Forecasting: Predicting the actions of the firms' investors
- Methods
    - A simple logistic model from 1994
    - A better logistic model from 2012
    - A LASSO model including firms' disclosure text
    - [Optional] e**X**treme **G**radient **Boost**ing (XGBoost)


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

- `package:coefplot`
- `package:glmnet`
- `package:kableExtra`
- `package:knitr`
- `package:magrittr`
- `package:revealjs`
- `package:tidyverse`
- `package:yardstick`

## References {data-background="../Backgrounds/books-old.jpg" class="default stack never"}

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

<small style="width:100%">
```{r, eval=F}

```
</small>


# Appendix on `package:parsnip` with LASSO {data-background="../Backgrounds/ending.jpg" class="header"}

## LASSO using `tidymodels` {data-background="../Backgrounds/default.jpg" class="default"}

- There are many convenience packages in R to simplify workflows
    - `tidymodels` is a collection of such packages
        - `package:parsnip` helps run models on many different backends
        - `package:recipes` helps process and prep data 
        - `package:rsample` for cross validation
        - `package:workflows` to tie it all together

> We will use `tidymodels` to run a LASSO and an XGBoost model for misreporting detection

- Jared Lander gave a good talk on using tidy models, <a href="https://jaredlander.com/content/2018/11/ManyWaysToLasso2.html">Many ways To Lasso</a>, at DSSG

## Data prep with `package:recipes` {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, include=F}
BCE_eq <- as.formula(paste("AAER ~ logtotasset + rsst_acc + chg_recv + chg_inv +
  soft_assets + pct_chg_cashsales + chg_roa + issuance +
  oplease_dum + book_mkt + lag_sdvol + merger + bigNaudit +
  midNaudit + cffin + exfin + restruct + bullets + headerlen +
  newlines + alltags + processedsize + sentlen_u + wordlen_s +
  paralen_s + repetitious_p + sentlen_s + typetoken +
  clindex + fog + active_p + passive_p + lm_negative_p +
  lm_positive_p + allcaps + exclamationpoints + questionmarks + ",
  paste(paste0("Topic_",1:30,"_n_oI"), collapse=" + "), collapse=""))

```

```{r, message=F}
library(recipes)
library(parsnip)

df <- read_csv("../../Data/Session_6.csv")
BCEformula <- BCE_eq

train <- df %>% filter(Test == 0)
test <- df %>% filter(Test == 1)

rec <- recipe(BCEformula, data = train) %>%
  step_zv(all_predictors()) %>%  # Drop any variables with zero variance
  step_center(all_predictors()) %>%  # Center all prediction variables
  step_scale(all_predictors()) %>%  # Scale all prediction variables
  step_intercept() %>%  # Add an intercept to the model
  step_num2factor(all_outcomes(), ordered = T, levels=c("0","1"),
                  transform = function(x) x + 1)  # Convert DV to factor

prepped <- rec %>% prep(training=train)
```

## Running a model with `package:parsnip` {data-background="../Backgrounds/coding.jpg" class="default"}

```{r}
# "bake" your recipe to get data ready
train_baked  <- bake(prepped, new_data = train)
test_baked  <- bake(prepped, new_data = test)

# Run the model with parsnip
train_model <- logistic_reg(mixture=1, penalty=1) %>%  # mixture = 1 sets LASSO
  set_engine('glmnet') %>%
  fit(BCEformula, data = train_baked)
```

## Visualizing `package:parsnip`'s output {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, warning=F, message=F, fig.height=5}
# train_model$fit is the same as fit_LASSO earlier in the slides
coefplot(train_model$fit, lambda=0.002031, sort='magnitude')
```

## Plugging in to cross validation {data-background="../Backgrounds/coding.jpg" class="default"}

- `package:parsnip` can plug into cross validation through `package:rsample`, usingthrough `vfold_cv()`
    - Easy to do surface level analysis with it
    - Difficult to do anything more in depth still
- We can `juice()` out our data and just use `cv.glmnet()`

```{r}
rec <- recipe(BCEformula, data = train) %>%
  step_zv(all_predictors()) %>%  # Drop any variables with zero variance
  step_center(all_predictors()) %>%  # Center all prediction variables
  step_scale(all_predictors()) %>%  # Scale all prediction variables
  step_intercept()  # Add an intercept to the model

prepped <- rec %>% prep(training=train)
test_prepped <- rec %>% prep(training=test)

# "Juice" your recipe to get data for other packages
train_x <- juice(prepped, all_predictors(), composition = "dgCMatrix")
train_y <- juice(prepped, all_outcomes(), composition = "matrix")
test_x <- juice(test_prepped, all_predictors(), composition = "dgCMatrix")
test_y <- juice(test_prepped, all_outcomes(), composition = "matrix")
```

## Running a cross validated model {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, message=F}
# Cross validation
set.seed(75347)  #for reproducibility
cvfit = cv.glmnet(x=train_x, y=train_y, family = "binomial", alpha = 1,
                  type.measure="auc")
```

<div class="clearfix">
<div class="left">
```{r}
plot(cvfit)
```
</div>
<div class="right">
```{r}
cvfit$lambda.min
cvfit$lambda.1se
```
</div>
</div>

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

<div class="clearfix">
<div class="left">
`lambda.min`
```{r, warning=F, fig.height=5.2, fig.width=4, echo=F}
#coef(cvfit, s = "lambda.min")
coefplot(cvfit, lambda='lambda.min', sort='magnitude') + theme(axis.text.y = element_text(size=15))
```
</div>
<div class="right">
`lambda.1se`
```{r, warning=F, fig.height=5.2, fig.width=4, echo=F}
#coef(cvfit, s = "lambda.1se")
coefplot(cvfit, lambda='lambda.1se', sort='magnitude') + theme(axis.text.y = element_text(size=15))
```
</div>
</div>

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

```{r, echo=F, fig.height=3.5}
train.pred_L1.min2 <- c(predict(cvfit, train_x, type="response", s = "lambda.min"))
test.pred_L1.min2 <- c(predict(cvfit, test_x, type="response", s = "lambda.min"))
train.pred_L1.1se2 <- c(predict(cvfit, train_x, type="response", s = "lambda.1se"))
test.pred_L1.1se2 <- c(predict(cvfit, test_x, type="response", s = "lambda.1se"))

df_train <- data.frame(train.pred_L1.min2=train.pred_L1.min2, train.pred_L1.1se2=train.pred_L1.1se2, AAER=factor(train_y, levels=c(0,1)))
df_test <- data.frame(test.pred_L1.min2=test.pred_L1.min2, test.pred_L1.1se2=test.pred_L1.1se2, AAER=factor(test_y, levels=c(0,1)))

auc_in <- df_train %>% roc_auc(truth=AAER, estimate=train.pred_L1.min2, event_level='second')
auc_out <- df_test %>% roc_auc(truth=AAER, estimate=test.pred_L1.min2, event_level='second')
curve_in_CV.min2 <- df_train %>% roc_curve(truth=AAER, estimate=train.pred_L1.min2, event_level='second')
curve_out_CV.min2 <- df_test %>% roc_curve(truth=AAER, estimate=test.pred_L1.min2, event_level='second')

aucs_CV.min2 <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_CV.min2) <- c("In sample, lambda.min", "Out of sample, lambda.min")

auc_in <- df_train %>% roc_auc(truth=AAER, estimate=train.pred_L1.1se2, event_level='second')
auc_out <- df_test %>% roc_auc(truth=AAER, estimate=test.pred_L1.1se2, event_level='second')
curve_in_CV.1se2 <- df_train %>% roc_curve(truth=AAER, estimate=train.pred_L1.1se2, event_level='second')
curve_out_CV.1se2 <- df_test %>% roc_curve(truth=AAER, estimate=test.pred_L1.1se2, event_level='second')

aucs_CV.1se2 <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_CV.1se2) <- c("In sample, lambda.1se", "Out of sample, lambda.1se")

ggplot() +
  geom_line(data=curve_in_CV.min2, aes(y=sensitivity, x=1-specificity, color="In Sample, lambda.min")) + 
  geom_line(data=curve_out_CV.min2, aes(y=sensitivity, x=1-specificity, color="Out of Sample, lambda.min")) +
  geom_line(data=curve_in_CV.1se2, aes(y=sensitivity, x=1-specificity, color="In Sample, lambda.1se")) + 
  geom_line(data=curve_out_CV.1se2, aes(y=sensitivity, x=1-specificity, color="Out of Sample, lambda.1se")) +
  geom_abline(slope=1)

c(aucs_CV.min2, aucs_CV.1se2)
```

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

- `package:glmnet`
- `package:parsnip`
- `package:recipes`
- `package:yardstick`

# If you really want to use `package:parsnip` for CV LASSO {data-background="../Backgrounds/ending.jpg" class="header"}

## Data prep with `package:recipes` (Same as before) {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, include=F}
BCE_eq <- as.formula(paste("AAER ~ logtotasset + rsst_acc + chg_recv + chg_inv +
  soft_assets + pct_chg_cashsales + chg_roa + issuance +
  oplease_dum + book_mkt + lag_sdvol + merger + bigNaudit +
  midNaudit + cffin + exfin + restruct + bullets + headerlen +
  newlines + alltags + processedsize + sentlen_u + wordlen_s +
  paralen_s + repetitious_p + sentlen_s + typetoken +
  clindex + fog + active_p + passive_p + lm_negative_p +
  lm_positive_p + allcaps + exclamationpoints + questionmarks + ",
  paste(paste0("Topic_",1:30,"_n_oI"), collapse=" + "), collapse=""))

```

```{r, message=F}
library(tidyr)
library(tidymodels)
library(tidyverse)

df <- read_csv("../../Data/Session_6.csv")
BCEformula <- BCE_eq

train <- df %>% filter(Test == 0)
test <- df %>% filter(Test == 1)

LASSO_rec <- recipe(BCEformula, data = train) %>%
  step_zv(all_predictors()) %>%  # Drop any variables with zero variance
  step_center(all_predictors()) %>%  # Center all prediction variables
  step_scale(all_predictors()) %>%  # Scale all prediction variables
  step_intercept() %>%  # Add an intercept to the model
  step_num2factor(all_outcomes(), ordered = T, levels=c("0","1"),
                  transform = function(x) x + 1)  # Convert DV to factor
```

## Define a tuning with `package:tune` and `package:tidyr` {data-background="../Backgrounds/coding.jpg" class="default"}

```{r}
LASSO_mod <- logistic_reg(penalty=tune(), mixture=1) %>%  # mixture = 1 sets LASSO
  set_engine('glmnet')

# Define a grid to tune over
grid <- expand_grid(penalty = exp(seq(-11,-4, length.out=100)))
```

- `tune()` replaces any parameters you would like to tune over
- Unlike with `cv.glmnet()`, we'll need to specify the range to tune over
    - The `expand_grid()` function from `package:tidyr` makes this easy
    - The `exp(seq())` part is to emulate `cv.glmnet()`'s tuning behavior

## Define a workflow with `package:workflows` {data-background="../Backgrounds/coding.jpg" class="default"}

```{r}
LASSO_wfl <- workflow() %>%
  add_model(LASSO_mod) %>%
  add_recipe(LASSO_rec)
```

> A workflow tells the various fitting and tuning functions in `package:tune` how to handle the data.  In other words, this will combine our model and recipe into 1 object.

## Run the model using `package:rsample`, `package:tune`, and `package:yardstick` {data-background="../Backgrounds/coding.jpg" class="default"}

```{r}
set.seed(354351)
folds <- vfold_cv(train, v=10)  # from rsample
metrics = metric_set(roc_auc)  # from yardstick

LASSO_fit_tuned <- tune_grid(LASSO_wfl,
                             grid = grid,
                             resamples = folds,
                             metrics=metrics)
```

## Take a look at the output {data-background="../Backgrounds/coding.jpg" class="default"}

```{r}
LASSO_fit_tuned %>%
  collect_metrics()
```

## Plotting it out {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, fig.height=3.0}
lambda.min <- LASSO_fit_tuned %>%
  collect_metrics() %>%
  arrange(-mean) %>%
  slice(1) %>%
  pull(penalty) %>%
  log()

LASSO_fit_tuned %>%
  collect_metrics() %>%
  ggplot(aes(x=log(penalty), y=mean)) +
  geom_point() +
  xlab("Log(lambda)") + 
  geom_vline(xintercept = lambda.min)
```

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

- `package:glmnet`
- `package:parsnip`
- `package:recipes`
- `package:rsample`
- `package:tidyr`
- `package:tune`
- `package:workflows`
- `package:yardstick`

# Appendix on XGBoost {data-background="../Backgrounds/ending.jpg" class="header"}

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

- e**X**treme **G**radient **B**oosting
- A simple explanation:
    1. Start with 1 or more decision trees & check error
    2. Make more decision trees & check error
    3. Use the difference in error to guess a another model
    4. Repeat #2 and #3 until the model's error is stable

<img src="../Figures/XGBoost.png" width="600px">

## Data prep with `package:recipes` {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, include=F}
BCE_eq <- as.formula(paste("AAER ~ logtotasset + rsst_acc + chg_recv + chg_inv +
  soft_assets + pct_chg_cashsales + chg_roa + issuance +
  oplease_dum + book_mkt + lag_sdvol + merger + bigNaudit +
  midNaudit + cffin + exfin + restruct + bullets + headerlen +
  newlines + alltags + processedsize + sentlen_u + wordlen_s +
  paralen_s + repetitious_p + sentlen_s + typetoken +
  clindex + fog + active_p + passive_p + lm_negative_p +
  lm_positive_p + allcaps + exclamationpoints + questionmarks + ",
  paste(paste0("Topic_",1:30,"_n_oI"), collapse=" + "), collapse=""))

```

```{r, message=F}
library(recipes)
library(parsnip)

df <- read_csv("../../Data/Session_6.csv")
BCEformula <- BCE_eq

train <- df %>% filter(Test == 0)
test <- df %>% filter(Test == 1)

rec <- recipe(BCEformula, data = train) %>%
  step_zv(all_predictors()) %>%  # Drop any variables with zero variance
  step_center(all_predictors()) %>%  # Center all prediction variables
  step_scale(all_predictors()) %>%  # Scale all prediction variables
  step_intercept()  # Add an intercept to the model
```

```{r, message=F}
# Juice our data
prepped <- rec %>% prep(training=train)
train_x <- juice(prepped, all_predictors(), composition = "dgCMatrix")
train_y <- juice(prepped, all_outcomes(), composition = "matrix")
test_prepped <- rec %>% prep(training=test)
test_x <- juice(test_prepped, all_predictors(), composition = "dgCMatrix")
test_y <- juice(test_prepped, all_outcomes(), composition = "matrix")
```

## Running a cross validated model {data-background="../Backgrounds/coding.jpg" class="default"}


<div class="clearfix">
<div class="left">
```{r, message=F, warning=F, error=F}
# Cross validation
set.seed(482342)  #for reproducibility
library(xgboost)

# model setup
params <- list(max_depth=10,
               eta=0.2,
               gamma=10,
               min_child_weight = 5,
               objective =
                 "binary:logistic")

# run the model
xgbCV <- xgb.cv(params=params,
                data=train_x,
                label=train_y,
                nrounds=100,
                eval_metric="auc",
                nfold=10,
                stratified=TRUE)
```
</div>
<div class="right">
```{r}
numTrees <- min(
 which(
  xgbCV$evaluation_log$test_auc_mean == 
  max(xgbCV$evaluation_log$test_auc_mean)
 )
)

fit4 <- xgboost(params=params,
                data = train_x,
                label = train_y,
                nrounds = numTrees,
                eval_metric="auc")
```
</div>
</div>

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

```{r, fig.height=4.5}
xgb.train.data = xgb.DMatrix(train_x, label = train_y, missing = NA)
col_names = attr(xgb.train.data, ".Dimnames")[[2]]
imp = xgb.importance(col_names, fit4)
# Variable importance
xgb.plot.importance(imp)
```


## Model comparison: Out of sample {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, echo=F, fig.height=3.5}
pred.xgb <- c(predict(fit4, test_x, type="response"))
df_test <- data.frame(pred.xgb=pred.xgb, AAER=factor(test_y, levels=c(0,1)))

pred.xgb <- c(predict(fit4, train_x, type="response"))
df_train <- data.frame(pred.xgb=pred.xgb, AAER=factor(train_y, levels=c(0,1)))

auc_in <- df_train %>% roc_auc(truth=AAER, estimate=pred.xgb, event_level='second')
auc_out <- df_test %>% roc_auc(truth=AAER, estimate=pred.xgb, event_level='second')
curve_in_xgb <- df_train %>% roc_curve(truth=AAER, estimate=pred.xgb, event_level='second')
curve_out_xgb <- df_test %>% roc_curve(truth=AAER, estimate=pred.xgb, event_level='second')

aucs_xgb <- c(auc_in$.estimate, auc_out$.estimate)
names(aucs_xgb) <- c("In sample, XGBoost", "Out of sample, XGBoost")

ggplot() +
  geom_line(data=curve_out_BCE, aes(y=sensitivity, x=1-specificity, color="BCE")) + 
  geom_line(data=curve_out_2000f, aes(y=sensitivity, x=1-specificity, color="2000s + 2011")) + 
  geom_line(data=curve_out_2000s, aes(y=sensitivity, x=1-specificity, color="2000s")) + 
  geom_line(data=curve_out_2011, aes(y=sensitivity, x=1-specificity, color="2011")) + 
  geom_line(data=curve_out_1990s, aes(y=sensitivity, x=1-specificity, color="1990s")) +
  geom_line(data=curve_out_CV.min, aes(y=sensitivity, x=1-specificity, color="LASSO, lambda.min")) + 
  geom_line(data=curve_out_CV.1se, aes(y=sensitivity, x=1-specificity, color="LASSO, lambda.1se")) + 
  geom_line(data=curve_out_xgb, aes(y=sensitivity, x=1-specificity, color="XGBoost")) +
  geom_abline(slope=1)

oos_aucs <- c(aucs_1990s[2], aucs_2000s[2], aucs_2000f[2], aucs_2011[2], aucs_BCE[2], aucs_CV.1se[2], aucs_CV.min[2], auc_out$.estimate)
names(oos_aucs) <- c('1990s', '2000s', '2000s + 2011', '2011', 'BC', 'LASSO, lambda.1se', 'LASSO, lambda.min', 'XGBoost')
oos_aucs
```


```{r, include=F, echo=F, warning=F, message=F}
df_save_0 <- df_save %>% filter(Test==0)
df_save_1 <- df_save %>% filter(Test==1)
df_save_0$pred_xgb <- predict(fit4, train_x, type="response")
df_save_1$pred_xgb <- predict(fit4, test_x, type="response")
df_save <- bind_rows(df_save_0, df_save_1)
saveRDS(df_save, '../../Data/Session_6_models.rds')
```

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

- `package:parsnip`
- `package:recipes`
- `package:xgboost`
- `package:yardstick`

## NOTES {class="never"}

- Censoring is still difficult for students to wrap their head around