---
title: "Exercise template for 'Prediction'"
author: Your name
format: html
bibliography: book.bib
---

```{r}
#| message: false
#| include: false
library(tidyverse)
library(DBI)
library(farr)
library(furrr)
library(rpart)
library(rpart.plot)
library(glmnet)
```

```{r}
#| label: features
#| cache: false
#| include: false
db <- dbConnect(RPostgres::Postgres(), bigint = "integer")

funda <- tbl(db, Id(schema = "comp", table = "funda"))

X_vars <- c("act", "ap", "at", "ceq", "che", "cogs", "csho", "dlc", 
            "dltis", "dltt", "dp", "ib", "invt", "ivao", "ivst", 
            "lct", "lt", "ni", "ppegt", "pstk", "re", "rect", "sale", 
            "sstk", "txp", "txt", "xint", "prcc_f")

y_var <- "misstate"

features_all <-
  funda |>
  filter(indfmt == "INDL", datafmt == "STD",
         consol == "C", popsrc == "D") |>
  filter(fyear >= 1991, fyear <= 2008) |>
  mutate(across(c(ivao, pstk, ivst, txp), 
                \(x) coalesce(x, 0))) |>
  select(gvkey, datadate, fyear, all_of(X_vars))

features <-
  features_all |>
  arrange(gvkey, datadate) |>
  collect() |>
  na.omit()

dbDisconnect(db)
```

```{r}
#| include: false
aaer_long <-
  aaer_firm_year |>
  rowwise() |>
  mutate(fyear = list(seq(min_year, max_year, by = 1))) |>
  unnest(fyear) |>
  select(gvkey, p_aaer, fyear, min_year, max_year) 

get_aaers <- function(test_date, include_latent = FALSE) {
  df <-
    aaer_dates |>
    mutate(p_aaer = str_replace(aaer_num, "^AAER-", "")) |>
    select(p_aaer, aaer_date) |>
    distinct() |>
    right_join(aaer_long, by = "p_aaer") |>
    group_by(gvkey, fyear) |>
    summarize(misstate = any(aaer_date <= as.Date(test_date), na.rm = TRUE),
              .groups = "drop")

  if (include_latent) {
    df |> mutate(misstate_latent = TRUE)
  } else {
    df
  }
}
```

```{r}
#| label: prediction-models
#| cache: true
#| cache.vars: [data_train, n_folds, folds, sample_splits, formula, fm1, within_sample1]
#| cache.extra: !expr digest::digest(digest::digest(list(features, aaer_firm_year, aaer_dates, body(get_aaers))))
#| include: false
#| warning: false
aaers <- get_aaers(test_date = "2002-09-30", include_latent = TRUE)

df <-
  aaers |>
  select(-misstate_latent) |>
  right_join(features, by = c("gvkey", "fyear")) |>
  mutate(misstate = coalesce(misstate, FALSE)) 

data_train <- 
  df |>
  filter(fyear <= 2001)

n_folds <- 5
folds <- 1:n_folds

set.seed(2024)

sample_splits <- 
  data_train |>
  select(gvkey) |>
  distinct() |>
  arrange(gvkey) |>
  mutate(fold = sample(folds, nrow(pick(everything())), replace = TRUE))

formula <- str_c(y_var, " ~ ", str_c(X_vars, collapse = " + "))
fm1 <- glm(formula, data = data_train, family = binomial)

within_sample1 <-
  data_train |>
  mutate(score = predict(fm1, 
                         newdata = pick(everything()), type = "response"),
         predicted = as.integer(score > 0.5))
```

```{r}
#| include: false
logit_predict <- function(fold) {
  dft <-
    data_train |>
    inner_join(sample_splits, by = "gvkey")
  
  fm1 <-
    dft |>
    filter(fold != !!fold) |>
    glm(formula, data = _, family = binomial)
  
  dft |>
    filter(fold == !!fold) |>
    mutate(score = predict(fm1, pick(everything()), type = "response"),
           predicted = as.integer(score > 0.5)) |>
    select(gvkey, fyear, score, predicted)
}
```

```{r}
#| include: false
#| warning: false
#| cache: true
#| label: logit-cv
#| dependson: prediction-models
#| cache.vars: [logit_fit]
#| cache.extra: !expr digest::digest(body(logit_predict))
logit_fit <- 
  folds |> 
  map(logit_predict) |>
  list_rbind() |>
  inner_join(data_train, by = c("gvkey", "fyear")) |> 
  select(gvkey, fyear, predicted, score, misstate)
```

```{r}
#| label: tbl-logit-baseline
#| tbl-cap: "Supplied logistic regression performance"
#| echo: false
list(`In sample` = within_sample1, `Firm-fold CV` = logit_fit) |>
  imap_dfr(\(d, label) tibble(
    sample = label, n = nrow(d),
    auc = unname(auc(d$score, d$misstate)),
    true_positive = sum(d$predicted == 1 & d$misstate),
    false_positive = sum(d$predicted == 1 & !d$misstate),
    false_negative = sum(d$predicted == 0 & d$misstate))) |>
  knitr::kable(digits = 4)
```


### Classification trees

```{r}
#| include: false
fm2 <- rpart(formula, data = data_train, 
             control = rpart.control(cp = 0.001, minbucket = 5, xval = 0),
             method = "class")
```

```{r}
#| include: false
within_sample2 <-
  data_train |>
  mutate(score = predict(fm2, newdata = pick(everything()))[, "TRUE"],
         predicted = score > 0.5)
```

```{r}
#| echo: false
#| label: fig-tree
#| fig-height: 3
#| fig-cap: Fitted tree produced by recursive partitioning.
#| fig-alt: "Diagram of the fitted classification tree, showing split conditions and the class assigned at each terminal node."
rpart.plot(fm2, extra = 1)
```

### Exercises

1. One claimed benefit of classification trees is their ease of interpretation.
Can you provide an intuitive explanation for the fitted tree in @fig-tree?
If so, outline your explanation.
If not, what challenges prevent you from doing so?

2. Use `r n_folds`-fold cross-validation to assess the performance of the classification tree approach using the parameters above.

3. What do the parameters `cp` and `minbucket` represent?

4. How does the classification tree change if you increase the parameter `cp`?

5. Use `r n_folds`-fold cross-validation and AUC to assess the performance of the classification tree approach with three different parameters for `minbucket` (the one used above and two others).

6. Calculate AUC for the data stored in `within_sample2` above (this is for the tree stored in `fm2` and shown in @fig-tree).

```{r}
#| include: false
#| eval: false
auc(within_sample2$score, within_sample2$misstate)
```

## Penalized models {#sec-predict-penalized}

```{r}
#| include: false
#| label: fm-lasso
#| cache: true
#| cache.vars: [dft, fm_lasso_cv, idmin, fit_lasso_cv, fm_lasso]
#| dependson: prediction-models
#| cache.extra: !expr digest::digest(list(body(logit_predict), if (exists("rpart_predict")) body(rpart_predict), packageVersion("rpart"), packageVersion("glmnet")))
dft <-
  data_train |>
  inner_join(sample_splits, by = "gvkey")

fm_lasso_cv <-
  cv.glmnet(x = as.matrix(dft[X_vars]),
            y = dft[[y_var]],
            family = "binomial",
            alpha = 1,
            type.measure = "auc",
            foldid = dft[["fold"]],
            keep = TRUE)

idmin <- match(fm_lasso_cv$lambda.min, fm_lasso_cv$lambda)

fit_lasso_cv <- 
  dft |>
  select(misstate) |>
  mutate(logodds = fm_lasso_cv$fit.preval[, idmin],
         prob = exp(logodds)/(1 + exp(logodds)),
         predicted = logodds > 0)

fm_lasso <- glmnet(x = as.matrix(data_train[X_vars]), 
                   y = data_train[[y_var]],
                   family = "binomial",
                   lambda = fm_lasso_cv$lambda.min,
                   alpha = 1)
```

```{r}
#| label: tbl-lasso-baseline
#| tbl-cap: "Supplied lasso cross-validation performance"
#| echo: false
lasso_id <- match(fm_lasso_cv$lambda.min, fm_lasso_cv$lambda)
tibble(metric = c("Mean fold AUC", "Pooled out-of-fold AUC"),
       value = c(fm_lasso_cv$cvm[lasso_id],
                 unname(auc(fit_lasso_cv$prob, fit_lasso_cv$misstate)))) |>
  knitr::kable(digits = 4)
```


### Discussion questions and exercises

1. What features are selected by the lasso model in this case?
(*Hint*: You may find it useful to inspect `fm_lasso$beta[, 1]`.)

2. Modify the code above to estimate a ridge regression model.
(*Hint*: Set `alpha = 0`.)

3. Describe how would you estimate an elastic net model.
(*Hint*: The help page for `cv.glmnet()` from the `glmnet` package might be useful. One approach would be to put some of the steps above in a function that accepts `alpha` as an argument and then estimate over different values of $\alpha$.)

4. Calculate AUC and produce the confusion matrix for `fm_lasso` applied to `data_train` (i.e., in sample).
Interpret these results.
(*Hint*: You might use this code: `predict(fm_lasso, newx = as.matrix(data_train[X_vars]))` and the `auc()` function supplied with the `farr` package.)

```{r}
#| label: fit_rus_model
#| cache: false
#| include: false
fit_rus_model <- function(df, size = 30, rus = TRUE, learn_rate = 1,
                          maxdepth = NULL, minbucket = NULL,
                          ir = 1) {
  control_args <- list()
  if (!is.null(maxdepth)) control_args$maxdepth <- maxdepth
  if (!is.null(minbucket)) control_args$minbucket <- minbucket
  control <- do.call(rpart.control, control_args)
  
  fm <- rusboost(formula, df, size = size, ir = ir, learn_rate = learn_rate,
                 rus = rus, control = control)
  return(fm)
}
```

```{r}
#| include: false
#| label: rus_predict
#| cache: false
rus_predict <- function(fold, ...) {
  
  dft <-
    data_train |>
    inner_join(sample_splits, by = "gvkey") |>
    mutate(misstate = as.integer(misstate))
  
  fm <-
    dft |>
    filter(fold != !!fold) |>
    fit_rus_model(...)
  
  res <-
    dft |>
    filter(fold == !!fold) |>
    mutate(prob = predict(fm, pick(everything()), type = "prob"),
           predicted = predict(fm, pick(everything()), type = "class")) |>
    select(gvkey, fyear, prob, predicted)
  
  res
}
```

```{r}
#| include: false
#| label: get_auc
#| warning: false
#| cache: false
get_auc <- function(...) {
  set.seed(2021)
  rus_fit <- 
    folds |> 
    future_map(rus_predict, 
               .options = furrr_options(seed = 2021),
               ...) |>
    list_rbind() |>
    inner_join(data_train, by = c("gvkey", "fyear")) |> 
    select(gvkey, fyear, predicted, prob, misstate)
  
  auc(score = rus_fit$prob, response = as.integer(rus_fit$misstate))
}
```

```{r}
#| include: false
test_years <- 2003:2008

test_dates <-
  tibble(test_year = test_years) |>
  mutate(test_date = as.Date(paste0(test_year + 1, "-09-30")))
```

```{r}
#| include: false
train_model <- function(features, test_year, gap = 2,
                        size, rus = TRUE, learn_rate, ir, control) {

  test_date <-
    test_dates |>
    filter(test_year == !!test_year) |>
    select(test_date) |>
    pull()

  aaers <- get_aaers(test_date = test_date)

  data_train <-
    aaers |>
    right_join(features, by = c("gvkey", "fyear")) |>
    filter(fyear <= test_year - gap) |>
    mutate(misstate = as.integer(coalesce(misstate, FALSE)))

  fm <- rusboost(formula, data_train, size = size, learn_rate = learn_rate,
                 rus = rus, ir = ir, control = control)
  return(fm)
}
```

```{r}
#| include: false
test_model <- function(features, fm, test_year) {

  target <- as.character(stats::as.formula(fm[["formula"]])[[2]])

  df <-
    features |>
    left_join(aaer_long |> distinct(gvkey, fyear) |>
                mutate(misstate = 1L), by = c("gvkey", "fyear")) |>
    mutate(misstate = coalesce(misstate, 0L))

  data_test <-
    df |>
    filter(fyear == test_year)

  return(list(score = predict(fm, data_test, type = "prob"),
              predicted = as.integer(predict(fm, data_test,
                                             type = "class") == 1),
              response = as.integer(data_test[[target]] == 1)))
}
```

```{r}
#| include: false
fit_test <- function(x, features, gap = 2,
                     size, rus = TRUE, learn_rate, ir = 1, control) {
  set.seed(2021)
  fm <- train_model(features = features, test_year = x,
                    size = size, rus = rus, ir = ir,
                    learn_rate = learn_rate, control = control, gap = gap)
  test_model(features = features, fm = fm, test_year = x)
}
```

```{r}
#| include: false
get_stats <- function(results) {
  response <- results$response
  score <- results$score

  ndcg_k <- function(k) {
    c(ndcg(score, response, k = k))
  }

  res <- list(auc(score, response),
              unlist(map(seq(0.01, 0.03, 0.01), ndcg_k)))
  unlist(res)
}
```

```{r}
#| include: false
add_means <- function(df) {
  bind_rows(df, summarise(df, across(everything(), mean)))
}
```

```{r}
#| include: false
table_results <- function(results, ...) {
  results |>
    map(get_stats) |>
    bind_rows() |>
    add_means() |>
    mutate(test_year = c(test_years, "Mean")) |>
    select(test_year, everything())
}
```

```{r}
#| include: false
optim_ada <- tibble(size = 50, learn_rate = 0.1, maxdepth = 7)
optim_rus <- tibble(size = 300, learn_rate = 0.01, maxdepth = 15, ir = 2)
plan(multisession, workers = min(2L, future::availableCores()))
```

```{r}
#| label: supplied-boost-cv
#| include: false
#| cache: true
#| cache.vars: [boost_cv]
#| dependson: prediction-models
#| cache.extra: !expr digest::digest(list(optim_ada, optim_rus, body(get_auc), body(rus_predict), body(fit_rus_model)))
boost_cv <- tibble(
  algorithm = c("AdaBoost", "RUSBoost"),
  auc = c(get_auc(size = optim_ada$size, rus = FALSE,
                  learn_rate = optim_ada$learn_rate, maxdepth = optim_ada$maxdepth),
          get_auc(size = optim_rus$size, rus = TRUE, ir = optim_rus$ir,
                  learn_rate = optim_rus$learn_rate, maxdepth = optim_rus$maxdepth)))
```

```{r}
#| label: supplied-boost-tests
#| include: false
#| cache: true
#| cache.vars: [results_ada, results_rus]
#| dependson: prediction-models
#| cache.extra: !expr digest::digest(list(features, aaer_long, test_dates, optim_ada, optim_rus, body(train_model), body(test_model), body(fit_test)))
plan(multisession, workers = min(2L, future::availableCores()))
results_ada <- future_map(test_years, fit_test, features = features,
                          size = optim_ada$size, rus = FALSE,
                          learn_rate = optim_ada$learn_rate,
                          control = rpart.control(maxdepth = optim_ada$maxdepth),
                          .options = furrr_options(seed = 2021))
results_rus <- future_map(test_years, fit_test, features = features,
                          size = optim_rus$size, rus = TRUE, ir = optim_rus$ir,
                          learn_rate = optim_rus$learn_rate,
                          control = rpart.control(maxdepth = optim_rus$maxdepth),
                          .options = furrr_options(seed = 2021))
```

The supplied boosting settings are those embedded in the book's existing tuning output.
The cross-validation scores and chronological test scores below are recomputed on the current data using those settings; the full tuning grids are not rerun in these companions.

```{r}
#| label: tbl-boost-cv-supplied
#| tbl-cap: "Supplied boosting models: Current cross-validation AUC"
#| echo: false
boost_cv |> knitr::kable(digits = 4)
```

```{r}
#| label: tbl-ada-supplied
#| tbl-cap: "Supplied AdaBoost test-year performance"
#| echo: false
table_results(results_ada) |> knitr::kable(digits = 4)
```

```{r}
#| label: tbl-rus-supplied
#| tbl-cap: "Supplied RUSBoost test-year performance"
#| echo: false
table_results(results_rus) |> knitr::kable(digits = 4)
```

## Discussion questions and exercises

1. We did not evaluate the statistical significance of the difference between the AUCs for AdaBoost and RUSBoost.
Describe how you might evaluate the statistical significance of this difference.

2. Compare the performance statistics from cross-validation and out-of-sample testing.
What do these numbers tell you?
Do these provide evidence on whether using (say) 2001 data to train a model that is evaluated on (say) 1993 frauds is problematic?

3. @Bao:2022aa report an AUC in the test sample for RUSBoost of 0.7228 (see Panel B of Table 1).
What differences are there between our approach here and those in @Bao:2020aa and @Bao:2022aa that might account for the differences in AUC?

4. @Bao:2022aa report an AUC in the test sample for an all-features logit model of 0.6842 (see Panel B of Table 1).
Calculate an equivalent value for AUC from applying logistic regression to data from above.
You may find it helpful to adapt the code above for RUSBoost in the following steps:
   a. Create a function `train_logit()` analogous to `train_model()` above, but using `glm()` in place of `rusboost()`.
   b. Create a function `test_logit()` analogous to `test_model()` above, but using calls to `predict()` like those used above for logistic regression.
   c. Create a function `fit_test_logit()` analogous to `fit_test()` above, but calling the functions you created above.
   d. Use `map()` or `future_map()` with `test_years` and `fit_test_logit()` and store the results in `results_logit`.
   e. Tabulate results using `table_results(results_logit)`. ([the book's logit results](https://iangow.github.io/far_book/prediction.html#tbl-logit-oss) contains our results from following these steps.)
   
5. Do the AUC results obtained from logit (either from your answer to the previous question or [the book's logit results](https://iangow.github.io/far_book/prediction.html#tbl-logit-oss)) surprise you given the cross-validation and in-sample AUCs calculated (in exercises) above?
How does the average AUC compare with the 0.6842 of @Bao:2022aa?
What might account for the difference?

6. Provide an intuitive explanation of the out-of-sample (test) NDCG\@k results for the RUSBoost model.
(Imagine you are trying to explain these to a user of prediction models like those discussed earlier in this chapter.)
Do you believe that these results are strong, weak, or somewhere in between?

7. We followed @Bao:2020aa in using `gap = 2` (i.e., training data for fiscal years up to two years before the test year), but discussed the trade-offs in using a shorter gap period (more data, but more measurement error).
How would you evaluate alternative gap periods?

8. Which results are emphasized in @Bao:2020aa?
What are the equivalent values in @Bao:2022aa?
Which results are emphasized in @Bao:2022aa?

9. Which audience (e.g., practitioners, researchers) would be most interested in the results of @Bao:2020aa and @Bao:2022aa?
What limitations of @Bao:2020aa and @Bao:2022aa might affect the interest of this audience?

10. Provide reasons why the tree-based models outperform the logit-based models in this setting?
Do you think this might be related to the set of features considered?
How might you alter the features if you wanted to improve the performance of the logit-based models?

```{r}
#| include: false
test_years <- 2003:2008

test_dates <- 
  tibble(test_year = test_years) |>
  mutate(test_date = as.Date(paste0(test_year + 1, "-09-30")))
```