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

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

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

```{r, warning=FALSE}
library(tidyverse)
library(coefplot)
library(DT)
```

```{r}
df <- read.csv("../../Data/Session_6.csv")
```

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

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

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

```{r, warning=F, message=F, fig.height=5}
library(ROCR)
pred <- predict(fit_1990s, df, type="response")
ROCpred <- prediction(as.numeric(pred[df$Test==0 & !is.na(pred)]), as.numeric(df[df$Test==0 & !is.na(pred),]$AAER))
ROCpred_out <- prediction(as.numeric(pred[df$Test==1 & !is.na(pred)]), as.numeric(df[df$Test==1 & !is.na(pred),]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_1990s <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_1990s <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))

ggplot() +
  geom_line(data=df_ROC_1990s, aes(x=FalsePositive, y=TruePositive, color="In Sample")) +
  geom_line(data=df_ROC_out_1990s, aes(x=FalsePositive, y=TruePositive, color="Out of Sample")) + 
  geom_abline(slope=1)

auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_1990s <- c(auc@y.values[[1]], auc_out@y.values[[1]])
names(aucs_1990s) <- c("In sample AUC", "Out of sample AUC")
aucs_1990s
```

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

```{r, warning=F, message=F, fig.height=5}
library(ROCR)
pred <- predict(fit_2011, df, type="response")
ROCpred <- prediction(as.numeric(pred[df$Test==0 & !is.na(pred)]), as.numeric(df[df$Test==0 & !is.na(pred),]$AAER))
ROCpred_out <- prediction(as.numeric(pred[df$Test==1 & !is.na(pred)]), as.numeric(df[df$Test==1 & !is.na(pred),]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_2011 <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_2011 <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))

ggplot() +
  geom_line(data=df_ROC_2011, aes(x=FalsePositive, y=TruePositive, color="In Sample")) +
  geom_line(data=df_ROC_out_2011, aes(x=FalsePositive, y=TruePositive, color="Out of Sample")) + 
  geom_abline(slope=1)

auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_2011 <- c(auc@y.values[[1]], auc_out@y.values[[1]])
names(aucs_2011) <- c("In sample AUC", "Out of sample AUC")
aucs_2011
```

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

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

```{r, warning=F, message=F, fig.height=5}
library(ROCR)
pred <- predict(fit_2000s, df, type="response")
ROCpred <- prediction(as.numeric(pred[df$Test==0 & !is.na(pred)]), as.numeric(df[df$Test==0 & !is.na(pred),]$AAER))
ROCpred_out <- prediction(as.numeric(pred[df$Test==1 & !is.na(pred)]), as.numeric(df[df$Test==1 & !is.na(pred),]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_2000s <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_2000s <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))

ggplot() +
  geom_line(data=df_ROC_2000s, aes(x=FalsePositive, y=TruePositive, color="In Sample")) +
  geom_line(data=df_ROC_out_2000s, aes(x=FalsePositive, y=TruePositive, color="Out of Sample")) + 
  geom_abline(slope=1)

auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_2000s <- c(auc@y.values[[1]], auc_out@y.values[[1]])
names(aucs_2000s) <- c("In sample AUC", "Out of sample AUC")
aucs_2000s
```

```{r, warning=F, message=F}
df_save$pred_S <- pred
```

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

```{r, warning=F, message=F, fig.height=5}
library(ROCR)
pred <- predict(fit_2000f, df, type="response")
ROCpred <- prediction(as.numeric(pred[df$Test==0 & !is.na(pred)]), as.numeric(df[df$Test==0 & !is.na(pred),]$AAER))
ROCpred_out <- prediction(as.numeric(pred[df$Test==1 & !is.na(pred)]), as.numeric(df[df$Test==1 & !is.na(pred),]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_2000f <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_2000f <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))

ggplot() +
  geom_line(data=df_ROC_2000f, aes(x=FalsePositive, y=TruePositive, color="In Sample")) +
  geom_line(data=df_ROC_out_2000f, aes(x=FalsePositive, y=TruePositive, color="Out of Sample")) + 
  geom_abline(slope=1)

auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_2000f <- c(auc@y.values[[1]], auc_out@y.values[[1]])
names(aucs_2000f) <- c("In sample AUC", "Out of sample AUC")
aucs_2000f
```

```{r, warning=F, message=F}
df_save$pred_FS <- pred
```

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

```{r, warning=F, message=F, fig.height=5}
library(ROCR)
pred <- predict(fit_BCE, df, type="response")
ROCpred <- prediction(as.numeric(pred[df$Test==0 & !is.na(pred)]), as.numeric(df[df$Test==0 & !is.na(pred),]$AAER))
ROCpred_out <- prediction(as.numeric(pred[df$Test==1 & !is.na(pred)]), as.numeric(df[df$Test==1 & !is.na(pred),]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_BCE <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_BCE <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))

ggplot() +
  geom_line(data=df_ROC_BCE, aes(x=FalsePositive, y=TruePositive, color="In Sample")) +
  geom_line(data=df_ROC_out_BCE, aes(x=FalsePositive, y=TruePositive, color="Out of Sample")) + 
  geom_abline(slope=1)

auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_BCE <- c(auc@y.values[[1]], auc_out@y.values[[1]])
names(aucs_BCE) <- c("In sample AUC", "Out of sample AUC")
aucs_BCE
```

```{r, warning=F, message=F}
df_save$pred_BCE <- pred
```

```{r, warning=F, message=F, fig.height=5}
ggplot() +
  geom_line(data=df_ROC_out_BCE, aes(x=FalsePositive, y=TruePositive, color="BCE")) +
  geom_line(data=df_ROC_out_2000f, aes(x=FalsePositive, y=TruePositive, color="2000s + 2011")) + 
  geom_line(data=df_ROC_out_2000s, aes(x=FalsePositive, y=TruePositive, color="2000s")) + 
  geom_line(data=df_ROC_out_2011, aes(x=FalsePositive, y=TruePositive, color="2011")) + 
  geom_line(data=df_ROC_out_1990s, aes(x=FalsePositive, y=TruePositive, 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
```

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

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

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

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

```{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
pred <- predict(fit_LASSO, xp, type="response", s = 0.002031)
```

```{r, fig.height=4}

ROCpred <- prediction(as.numeric(pred[df$Test==0 & !is.na(pred)]), as.numeric(df[df$Test==0 & !is.na(pred),]$AAER))
ROCpred_out <- prediction(as.numeric(pred[df$Test==1 & !is.na(pred)]), as.numeric(df[df$Test==1 & !is.na(pred),]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_L1 <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_L1 <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))

ggplot() +
  geom_line(data=df_ROC_BCE, aes(x=FalsePositive, y=TruePositive, color="In Sample")) +
  geom_line(data=df_ROC_out_BCE, aes(x=FalsePositive, y=TruePositive, color="Out of Sample")) + 
  geom_abline(slope=1)

auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_L1 <- c(auc@y.values[[1]], auc_out@y.values[[1]])
names(aucs_L1) <- c("In sample AUC", "Out of sample AUC")
aucs_L1
```

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

```{r}
plot(cvfit)
```

```{r}
cvfit$lambda.min
cvfit$lambda.1se
```

```{r, warning=F, fig.height=5.2, fig.width=4}
#coef(cvfit, s = "lambda.min")
coefplot(cvfit, lambda='lambda.min', sort='magnitude') + theme(axis.text.y = element_text(size=15))
```

```{r, warning=F, fig.height=5.2, fig.width=4}
#coef(cvfit, s = "lambda.1se")
coefplot(cvfit, lambda='lambda.1se', sort='magnitude') + theme(axis.text.y = element_text(size=15))
```

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

```{r, warning=F, message=F}
df_save$pred_lmin <- pred
df_save$pred_l1se <- pred2
```

```{r, fig.height=3.5}
ROCpred <- prediction(as.numeric(pred[df$Test==0 & !is.na(pred)]), as.numeric(df[df$Test==0 & !is.na(pred),]$AAER))
ROCpred_out <- prediction(as.numeric(pred[df$Test==1 & !is.na(pred)]), as.numeric(df[df$Test==1 & !is.na(pred),]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_CV <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_CV <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))
auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_CV <- c(auc@y.values[[1]], auc_out@y.values[[1]])


ROCpred <- prediction(as.numeric(pred2[df$Test==0 & !is.na(pred)]), as.numeric(df[df$Test==0 & !is.na(pred),]$AAER))
ROCpred_out <- prediction(as.numeric(pred2[df$Test==1 & !is.na(pred)]), as.numeric(df[df$Test==1 & !is.na(pred),]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_CV2 <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_CV2 <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))
auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_CV2 <- c(auc@y.values[[1]], auc_out@y.values[[1]])


ggplot() +
  geom_line(data=df_ROC_CV, aes(x=FalsePositive, y=TruePositive, color="In Sample, lambda.min")) +
  geom_line(data=df_ROC_out_CV, aes(x=FalsePositive, y=TruePositive, color="Out of Sample, lambda.min")) +
  geom_line(data=df_ROC_CV2, aes(x=FalsePositive, y=TruePositive, color="In Sample, lambda.1se")) +
  geom_line(data=df_ROC_out_CV2, aes(x=FalsePositive, y=TruePositive, color="Out of Sample, lambda.1se")) + 
  geom_abline(slope=1)

aucs <- c(aucs_CV, aucs_CV2)
names(aucs) <- c("In sample AUC, lambda.min", "Out of sample AUC, lambda.min", "In sample AUC, lambda.1se", "Out of sample AUC, lambda.1se")
aucs
```

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

```{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) %>%  # mixture = 1 sets LASSO
  set_engine('glmnet') %>%
  fit(BCEformula, data = train_baked)
```

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

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

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

```{r}
plot(cvfit)
```

```{r}
cvfit$lambda.min
cvfit$lambda.1se
```

```{r, warning=F, fig.height=5.2, fig.width=4}
#coef(cvfit, s = "lambda.min")
coefplot(cvfit, lambda='lambda.min', sort='magnitude') + theme(axis.text.y = element_text(size=15))
```

```{r, warning=F, fig.height=5.2, fig.width=4}
#coef(cvfit, s = "lambda.1se")
coefplot(cvfit, lambda='lambda.1se', sort='magnitude') + theme(axis.text.y = element_text(size=15))
```

```{r, fig.height=3.5}
# s= specifies the version of the model to use
pred_train.min <- predict(cvfit, train_x, type="response", s = "lambda.min")
pred_train.1se <- predict(cvfit, train_x, type="response", s = "lambda.1se")
pred_test.min <- predict(cvfit, test_x, type="response", s = "lambda.min")
pred_test.1se <- predict(cvfit, test_x, type="response", s = "lambda.1se")

ROCpred <- prediction(as.numeric(pred_train.min), as.numeric(df[df$Test==0,]$AAER))
ROCpred_out <- prediction(as.numeric(pred_test.min), as.numeric(df[df$Test==1,]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_CV <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_CV <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))
auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_CV <- c(auc@y.values[[1]], auc_out@y.values[[1]])


ROCpred <- prediction(as.numeric(pred_train.1se), as.numeric(df[df$Test==0,]$AAER))
ROCpred_out <- prediction(as.numeric(pred_test.1se), as.numeric(df[df$Test==1,]$AAER))
ROCperf <- performance(ROCpred, 'tpr','fpr')
ROCperf_out <- performance(ROCpred_out, 'tpr','fpr')
df_ROC_CV2 <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                 TruePositive=c(ROCperf@y.values[[1]]))
df_ROC_out_CV2 <- data.frame(FalsePositive=c(ROCperf_out@x.values[[1]]),
                 TruePositive=c(ROCperf_out@y.values[[1]]))
auc <- performance(ROCpred, measure = "auc")
auc_out <- performance(ROCpred_out, measure = "auc")
aucs_CV2 <- c(auc@y.values[[1]], auc_out@y.values[[1]])


ggplot() +
  geom_line(data=df_ROC_CV, aes(x=FalsePositive, y=TruePositive, color="In Sample, lambda.min")) +
  geom_line(data=df_ROC_out_CV, aes(x=FalsePositive, y=TruePositive, color="Out of Sample, lambda.min")) +
  geom_line(data=df_ROC_CV2, aes(x=FalsePositive, y=TruePositive, color="In Sample, lambda.1se")) +
  geom_line(data=df_ROC_out_CV2, aes(x=FalsePositive, y=TruePositive, color="Out of Sample, lambda.1se")) + 
  geom_abline(slope=1)

aucs <- c(aucs_CV, aucs_CV2)
names(aucs) <- c("In sample AUC, lambda.min", "Out of sample AUC, lambda.min", "In sample AUC, lambda.1se", "Out of sample AUC, lambda.1se")
aucs

add_auc <- aucs_CV[2]
names(add_auc) <- c("LASSO, lambda.min")
oos_aucs <- c(oos_aucs, add_auc)
```

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

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

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

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

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

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

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

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

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

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

```{r, fig.height=3.5}
# Usual AUC calculation
pred.xgb <- predict(fit4, test_x, type="response")
ROCpred.xgb <- prediction(as.numeric(pred.xgb), as.numeric(test_y))
ROCperf.xgb <- performance(ROCpred.xgb, 'tpr','fpr')
#plot(ROCperf.xgb)
df_ROC.xgb <- data.frame(FalsePositive=c(ROCperf.xgb@x.values[[1]]),
                 TruePositive=c(ROCperf.xgb@y.values[[1]]))
ggplot() +
  geom_line(data=df_ROC_out_BCE, aes(x=FalsePositive, y=TruePositive, color="BCE")) +
  geom_line(data=df_ROC_out_2000f, aes(x=FalsePositive, y=TruePositive, color="2000s + 2011")) + 
  geom_line(data=df_ROC_out_2000s, aes(x=FalsePositive, y=TruePositive, color="2000s")) + 
  geom_line(data=df_ROC_out_2011, aes(x=FalsePositive, y=TruePositive, color="2011")) + 
  geom_line(data=df_ROC_out_1990s, aes(x=FalsePositive, y=TruePositive, color="1990s")) + 
  geom_line(data=df_ROC_out_CV, aes(x=FalsePositive, y=TruePositive, color="LASSO, lambda.min")) +
  geom_line(data=df_ROC.xgb, aes(x=FalsePositive, y=TruePositive, color="XGBoost")) + 
  geom_abline(slope=1) + ggtitle("ROC Curves across models")

auc.xgb <- performance(ROCpred.xgb, measure = "auc")
auc <- auc.xgb@y.values[[1]]
names(auc) <- c("XGBoost AUC")
oos_aucs <- c(oos_aucs, auc)
oos_aucs
```

```{r, 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 <- pred.xgb
df_save <- bind_rows(df_save_0, df_save_1)
saveRDS(df_save, '../../Data/Session_6_models.rds')
```

