---
title: "Session 5 R practice"
author: "Dr. Richard M. Crowley"
date: "ACCT 420, Session 5"
output:
  html_document:
    theme: spacelab
    highlight: default
    toc: true
    toc_float:
      collapsed: false
    self_contained: false
---

```{r, include=FALSE}
tutorial::go_interactive()
```

## Logistic models and ROC

To run your code, click *run*.  It will let you know if your code is correct or not, and will offer hints on what to fix if needed.  If you are stuck, you can see the solution by clicking *Solution*.

### Exercise 1: Running a logit regression

This is a continuation of Session_4's exercises.

```{r ex="create_b", type="pre-exercise-code"}
library(dplyr) # loads dplyr for data manipulation
df <- data.frame(tic=rep("WMT",194),
                 fyearq=c(1970,1970,1970,1970,1971,1971,1971,1971,1972,1972,1972,1972,1973,1973,1973,1973,1974,1974,1974,1974,1975,1975,1975,1975,1976,1976,1976,1976,1977,1977,1977,1977,1978,1978,1978,1978,1979,1979,1979,1979,1980,1980,1980,1980,1981,1981,1981,1981,1982,1982,1982,1982,1983,1983,1983,1983,1984,1984,1984,1984,1985,1985,1985,1985,1986,1986,1986,1986,1987,1987,1987,1987,1988,1988,1988,1988,1989,1989,1989,1989,1990,1990,1990,1990,1991,1991,1991,1991,1992,1992,1992,1992,1993,1993,1993,1993,1994,1994,1994,1994,1995,1995,1995,1995,1996,1996,1996,1996,1997,1997,1997,1997,1998,1998,1998,1998,1999,1999,1999,1999,2000,2000,2000,2000,2001,2001,2001,2001,2002,2002,2002,2002,2003,2003,2003,2003,2004,2004,2004,2004,2005,2005,2005,2005,2006,2006,2006,2006,2007,2007,2007,2007,2008,2008,2008,2008,2009,2009,2009,2009,2010,2010,2010,2010,2011,2011,2011,2011,2012,2012,2012,2012,2013,2013,2013,2013,2014,2014,2014,2014,2015,2015,2015,2015,2016,2016,2016,2016,2017,2017,2017,2017,2018,2018),
                 fqtr=c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2),
                 revtq=c(7.267,9.939,11.47,15.61,12.758,18.207,19.205,27.844,21.011,29.918,32.683,41.277,32.737,37.655,42.098,55.071,46.355,56.177,59.462,74.215,64.65,78.852,86.108,110.721,93.971,115.186,121.875,147.775,122.655,152.381,172.637,230.783,167.274,204.085,221.982,306.957,238.569,291.685,300.747,417.175,314.575,372.598,401.842,554.184,434.251,533.653,634.476,842.617,665.229,788.705,827.767,1094.551,854.734,1098.942,1166.594,1546.639,1234.768,1508.533,1583.573,2073.983,1655.661,1934.358,2087.533,2773.936,2344.253,2770.276,2949.008,3845.557,3211.761,3732.237,4015.017,5000.238,4282.012,4845.133,4991.113,6530.734,5373.258,6046.41,6283.496,8107.484,6768.191,7543.508,7930.949,10358.898,9280.566,10339.898,10627.5,13638.797,11649.398,13028.398,13683.797,17122,13920.398,16236.5,16826.898,20360.5,17686.098,19942.297,20417.699,24447.797,20440,22723,22914,27550,22772,25587,25644,30856,25409,28386,28777,35386,29819,33521,33509,40785,35129,38913,40899,51868,43447,46588,46181,57079,48565,53187,53185,64735,52126,56781,55765,66905,57224,63231,63036,75190,65434,70516,69282,82900,70755,76697,75397,89252,79675,85430,84467,99078,86378,92827,91794,107289,94940,102342,98345,108747,94214,100876,99373,113622,99811,103726,101952,116360,104189,109366,110226,123169,113010,114065,113800,127776,114070,116830,115688,129706,114960,120125,119001,131565,114826,120229,117408,129667,115904,120319,118179,130936,117542,122968,123179,136267,122690,128028),
                 recession=c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
df <- df %>% mutate(revtq_yoy = (revtq - lag(revtq, 4)) / lag(revtq, 4),
                    recession_lag = lag(recession),
                    revtq_down = ifelse(revtq_yoy<0,1,0),
                    revtq_yoy_lag = lag(revtq_yoy),
                    date = fyearq + (fqtr-1)/4)

```

```{r ex="create_b", type="sample-code"}
# The dataframe below contains historical revenue for Walmart since 1970
# All data is quarterly
# The file uses Compustat's nomenclature:
#     - revenue is revtq
#     - ticker is tic
#     - fyearq is year
#     - fqtr is quarter
#     - recession is an indicator for US recessions
#     - revtq_yoy is year over year revenue growth
#     - revtq_down is an indicator equal to 1 if revtq_yoy is negative, else 0
#     - revtq_yoy_lag is the lag of revtq_yoy
df[1:5,]

# Calculate a logistic model of a decrease in year over year quarterly revenue
# on the previous quarter's year over year revenue growth, store it in mod1
mod1 <- glm( ~ , data=, family=)

# Print out a summary of the model using summary
summary(mod1)

#END
```

```{r ex="create_b", type="solution"}
# The dataframe below contains historical revenue for Walmart since 1970
# All data is quarterly
# The file uses Compustat's nomenclature:
#     - revenue is revtq
#     - ticker is tic
#     - fyearq is year
#     - fqtr is quarter
#     - recession is an indicator for US recessions
#     - revtq_yoy is year over year revenue growth
#     - revtq_down is an indicator equal to 1 if revtq_yoy is negative, else 0
#     - revtq_yoy_lag is the lag of revtq_yoy
df[1:5,]

# Calculate a logistic model of a decrease in year over year quarterly revenue
# on the previous quarter's year over year revenue growth, store it in mod1
mod1 <- glm(revtq_down ~ revtq_yoy_lag, data=df, family=binomial)

# Print out a summary of the model using summary
summary(mod1)

#END
```

```{r ex="create_b", type="sct"}
# Template based on https://www.rdocumentation.org/packages/testwhat/versions/4.1.1
# Check if something is explicitly typed

test_student_typed('revtq_down ~ revtq_yoy_lag',not_typed_msg='Make sure to regress `revtq_down` on `revtq_yoy_lag`.')
test_student_typed('data=df',not_typed_msg='Make sure to using the `df` data frame for the test.')
test_student_typed('family=binomial',not_typed_msg='Make sure to set family to be `binomial`, so as to run a logistic regression.')
test_expression_output("mod1", incorrect_msg="Your regression model isn't quite right.")
test_student_typed('summary(mod1)',not_typed_msg='Your model is correct, but did you remember to summarize it with `summary()`?')

# test_student_typed('x <- 2', not_typed_msg='')

# Check if function was used in input code
# test_function('c',incorrect_msg='')  

# Requires an object `x` to have the same value as the solution
# test_object("x",incorrect_msg = "",undefined_msg = "")  

# Requires an onject with the same value of `x` in the solution
# test_an_object("x",undefined_msg="")

# Checks if output of student's code contains given evaluated expression
# test_output_contains("x",incorrect_msg = "")

# Check if a vector of predefined objects are unchanged
# test_predefined_objects(c('x','y'),incorrect_msg="Don't onverwrite the predefined variables")

# Checks for a regex pattern in trhe output
# test_output_regex(pattern,fixed=F, times=1, incorrect_msg='')

# Can check an arbitrary expression across both solution and student code
#test_expression_output("typeof(company_name)", incorrect_msg="Did you store textual data in `company_name`?")

test_error()
success_msg("Awesome!")

# Other functions to note:
#     - test_or(a,b) -- checks if either test a or test b pass
#     - test_ggplot() -- can check if plots are correct
#     - test_function() -- can also check included parameters
#     - test_loop() -- checking for and while loops
#     - test_library_function('package', not_called_msg='',incorrect_msg='')
#     - test_if_else() -- checking if statements
#     - test_expression_error() -- can check if functions are properly defined
#     - test_operator('operator',), not_called_msg='',incorrect_msg='')
#     - test_function_definition() -- rigorously check defined function
#     - test_data_frame() -- check if dataframe [columns] are equivalent
#     - test_function_result, test_expression_result
```

### Exercise 2: Exploring model predictions

```{r ex="create_c", type="pre-exercise-code"}
library(dplyr) # loads dplyr for data manipulation
library(ggplot2)
df <- data.frame(tic=rep("WMT",194),
                 fyearq=c(1970,1970,1970,1970,1971,1971,1971,1971,1972,1972,1972,1972,1973,1973,1973,1973,1974,1974,1974,1974,1975,1975,1975,1975,1976,1976,1976,1976,1977,1977,1977,1977,1978,1978,1978,1978,1979,1979,1979,1979,1980,1980,1980,1980,1981,1981,1981,1981,1982,1982,1982,1982,1983,1983,1983,1983,1984,1984,1984,1984,1985,1985,1985,1985,1986,1986,1986,1986,1987,1987,1987,1987,1988,1988,1988,1988,1989,1989,1989,1989,1990,1990,1990,1990,1991,1991,1991,1991,1992,1992,1992,1992,1993,1993,1993,1993,1994,1994,1994,1994,1995,1995,1995,1995,1996,1996,1996,1996,1997,1997,1997,1997,1998,1998,1998,1998,1999,1999,1999,1999,2000,2000,2000,2000,2001,2001,2001,2001,2002,2002,2002,2002,2003,2003,2003,2003,2004,2004,2004,2004,2005,2005,2005,2005,2006,2006,2006,2006,2007,2007,2007,2007,2008,2008,2008,2008,2009,2009,2009,2009,2010,2010,2010,2010,2011,2011,2011,2011,2012,2012,2012,2012,2013,2013,2013,2013,2014,2014,2014,2014,2015,2015,2015,2015,2016,2016,2016,2016,2017,2017,2017,2017,2018,2018),
                 fqtr=c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2),
                 revtq=c(7.267,9.939,11.47,15.61,12.758,18.207,19.205,27.844,21.011,29.918,32.683,41.277,32.737,37.655,42.098,55.071,46.355,56.177,59.462,74.215,64.65,78.852,86.108,110.721,93.971,115.186,121.875,147.775,122.655,152.381,172.637,230.783,167.274,204.085,221.982,306.957,238.569,291.685,300.747,417.175,314.575,372.598,401.842,554.184,434.251,533.653,634.476,842.617,665.229,788.705,827.767,1094.551,854.734,1098.942,1166.594,1546.639,1234.768,1508.533,1583.573,2073.983,1655.661,1934.358,2087.533,2773.936,2344.253,2770.276,2949.008,3845.557,3211.761,3732.237,4015.017,5000.238,4282.012,4845.133,4991.113,6530.734,5373.258,6046.41,6283.496,8107.484,6768.191,7543.508,7930.949,10358.898,9280.566,10339.898,10627.5,13638.797,11649.398,13028.398,13683.797,17122,13920.398,16236.5,16826.898,20360.5,17686.098,19942.297,20417.699,24447.797,20440,22723,22914,27550,22772,25587,25644,30856,25409,28386,28777,35386,29819,33521,33509,40785,35129,38913,40899,51868,43447,46588,46181,57079,48565,53187,53185,64735,52126,56781,55765,66905,57224,63231,63036,75190,65434,70516,69282,82900,70755,76697,75397,89252,79675,85430,84467,99078,86378,92827,91794,107289,94940,102342,98345,108747,94214,100876,99373,113622,99811,103726,101952,116360,104189,109366,110226,123169,113010,114065,113800,127776,114070,116830,115688,129706,114960,120125,119001,131565,114826,120229,117408,129667,115904,120319,118179,130936,117542,122968,123179,136267,122690,128028),
                 recession=c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
df <- df %>% mutate(revtq_yoy = (revtq - lag(revtq, 4)) / lag(revtq, 4),
                    recession_lag = lag(recession),
                    revtq_down = ifelse(revtq_yoy<0,1,0),
                    revtq_yoy_lag = lag(revtq_yoy),
                    date = fyearq + (fqtr-1)/4)
mod1 <- glm(revtq_down ~ revtq_yoy_lag, data=df, family=binomial)
```

```{r ex="create_c", type="sample-code"}
# The dataframe below contains historical revenue for Walmart since 1970
# All data is quarterly
# The file uses Compustat's nomenclature:
#     - revtq_yoy is year over year revenue growth
#     - revtq_down is an indicator equal to 1 if revtq_yoy is negative, else 0
#     - revtq_yoy_lag is the lag of revtq_yoy
#     - date is a fractional date
# mod1 from the last exercise is also loaded

# Use predict() to calculate the probability of decreasing year over year
# revenue for Walmart based on mod1
df$pred <- predict(mod1, df, type="response")

# Plot out the actual decreases and predicted decrease probabilities by year
# ggplot2 is already loaded
ggplot(data=df) + 
  geom_point(aes(y = , x = , color = "Actual")) + 
  geom_point(aes(y = , x = , color = "Predicted")) + 
  ylab("Decrease in revenue")


#END
```

```{r ex="create_c", type="solution"}
# The dataframe below contains historical revenue for Walmart since 1970
# All data is quarterly
# The file uses Compustat's nomenclature:
#     - revtq_yoy is year over year revenue growth
#     - revtq_down is an indicator equal to 1 if revtq_yoy is negative, else 0
#     - revtq_yoy_lag is the lag of revtq_yoy
#     - date is a fractional date
# mod1 from the last exercise is also loaded

# Use predict() to calculate the probability of decreasing year over year
# revenue for Walmart based on mod1
df$pred <- predict(mod1, df, type="response")

# Plot out the actual decreases and predicted decrease probabilities by date
# ggplot2 is already loaded
ggplot(data=df) + 
  geom_point(aes(y = revtq_down, x = date, color = "Actual")) + 
  geom_point(aes(y = pred, x = date, color = "Predicted")) + 
  ylab("Decrease in revenue")


#END
```

```{r ex="create_c", type="sct"}
# Template based on https://www.rdocumentation.org/packages/testwhat/versions/4.1.1
# Check if something is explicitly typed

test_student_typed('x = date',not_typed_msg='Make sure to regress `revtq_down` on `revtq_yoy_lag`.')
test_student_typed('y = revtq_down',not_typed_msg='Make sure to set `y = revtq_down` within `aes()` for the actual decrease in revenue series.')
test_student_typed('y = pred',not_typed_msg='Make sure to set `y = pred` within `aes()` for the predicted decrease in revenue series.')
test_student_typed('type = "response"',not_typed_msg='Make sure to set `type = "response"` when using `predict()`.')
test_expression_output("df$pred", incorrect_msg="Your predictions in `df$pred` don't seem quite right.")

# test_student_typed('x <- 2', not_typed_msg='')

# Check if function was used in input code
# test_function('c',incorrect_msg='')  

# Requires an object `x` to have the same value as the solution
# test_object("x",incorrect_msg = "",undefined_msg = "")  

# Requires an onject with the same value of `x` in the solution
# test_an_object("x",undefined_msg="")

# Checks if output of student's code contains given evaluated expression
# test_output_contains("x",incorrect_msg = "")

# Check if a vector of predefined objects are unchanged
# test_predefined_objects(c('x','y'),incorrect_msg="Don't onverwrite the predefined variables")

# Checks for a regex pattern in trhe output
# test_output_regex(pattern,fixed=F, times=1, incorrect_msg='')

# Can check an arbitrary expression across both solution and student code
#test_expression_output("typeof(company_name)", incorrect_msg="Did you store textual data in `company_name`?")

test_error()
success_msg("Awesome!")

# Other functions to note:
#     - test_or(a,b) -- checks if either test a or test b pass
#     - test_ggplot() -- can check if plots are correct
#     - test_function() -- can also check included parameters
#     - test_loop() -- checking for and while loops
#     - test_library_function('package', not_called_msg='',incorrect_msg='')
#     - test_if_else() -- checking if statements
#     - test_expression_error() -- can check if functions are properly defined
#     - test_operator('operator',), not_called_msg='',incorrect_msg='')
#     - test_function_definition() -- rigorously check defined function
#     - test_data_frame() -- check if dataframe [columns] are equivalent
#     - test_function_result, test_expression_result
```

### Exercise 3: ROC curve and AUC

Note: This isn't using the RCOR library directly, but instead implements a subset of the functionality of RCOR 1.0-5.

```{r ex="create_d", type="pre-exercise-code"}
library(dplyr) # loads dplyr for data manipulation
library(ggplot2)
df <- data.frame(tic=rep("WMT",194),
                 fyearq=c(1970,1970,1970,1970,1971,1971,1971,1971,1972,1972,1972,1972,1973,1973,1973,1973,1974,1974,1974,1974,1975,1975,1975,1975,1976,1976,1976,1976,1977,1977,1977,1977,1978,1978,1978,1978,1979,1979,1979,1979,1980,1980,1980,1980,1981,1981,1981,1981,1982,1982,1982,1982,1983,1983,1983,1983,1984,1984,1984,1984,1985,1985,1985,1985,1986,1986,1986,1986,1987,1987,1987,1987,1988,1988,1988,1988,1989,1989,1989,1989,1990,1990,1990,1990,1991,1991,1991,1991,1992,1992,1992,1992,1993,1993,1993,1993,1994,1994,1994,1994,1995,1995,1995,1995,1996,1996,1996,1996,1997,1997,1997,1997,1998,1998,1998,1998,1999,1999,1999,1999,2000,2000,2000,2000,2001,2001,2001,2001,2002,2002,2002,2002,2003,2003,2003,2003,2004,2004,2004,2004,2005,2005,2005,2005,2006,2006,2006,2006,2007,2007,2007,2007,2008,2008,2008,2008,2009,2009,2009,2009,2010,2010,2010,2010,2011,2011,2011,2011,2012,2012,2012,2012,2013,2013,2013,2013,2014,2014,2014,2014,2015,2015,2015,2015,2016,2016,2016,2016,2017,2017,2017,2017,2018,2018),
                 fqtr=c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2),
                 revtq=c(7.267,9.939,11.47,15.61,12.758,18.207,19.205,27.844,21.011,29.918,32.683,41.277,32.737,37.655,42.098,55.071,46.355,56.177,59.462,74.215,64.65,78.852,86.108,110.721,93.971,115.186,121.875,147.775,122.655,152.381,172.637,230.783,167.274,204.085,221.982,306.957,238.569,291.685,300.747,417.175,314.575,372.598,401.842,554.184,434.251,533.653,634.476,842.617,665.229,788.705,827.767,1094.551,854.734,1098.942,1166.594,1546.639,1234.768,1508.533,1583.573,2073.983,1655.661,1934.358,2087.533,2773.936,2344.253,2770.276,2949.008,3845.557,3211.761,3732.237,4015.017,5000.238,4282.012,4845.133,4991.113,6530.734,5373.258,6046.41,6283.496,8107.484,6768.191,7543.508,7930.949,10358.898,9280.566,10339.898,10627.5,13638.797,11649.398,13028.398,13683.797,17122,13920.398,16236.5,16826.898,20360.5,17686.098,19942.297,20417.699,24447.797,20440,22723,22914,27550,22772,25587,25644,30856,25409,28386,28777,35386,29819,33521,33509,40785,35129,38913,40899,51868,43447,46588,46181,57079,48565,53187,53185,64735,52126,56781,55765,66905,57224,63231,63036,75190,65434,70516,69282,82900,70755,76697,75397,89252,79675,85430,84467,99078,86378,92827,91794,107289,94940,102342,98345,108747,94214,100876,99373,113622,99811,103726,101952,116360,104189,109366,110226,123169,113010,114065,113800,127776,114070,116830,115688,129706,114960,120125,119001,131565,114826,120229,117408,129667,115904,120319,118179,130936,117542,122968,123179,136267,122690,128028),
                 recession=c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
df <- df %>% mutate(revtq_yoy = (revtq - lag(revtq, 4)) / lag(revtq, 4),
                    recession_lag = lag(recession),
                    revtq_down = ifelse(revtq_yoy<0,1,0),
                    revtq_yoy_lag = lag(revtq_yoy),
                    date = fyearq + (fqtr-1)/4)
mod1 <- glm(revtq_down ~ revtq_yoy_lag, data=df, family=binomial)
df$pred <- predict(mod1, df, type="response")


prediction <- function(predictions, labels, label.ordering=NULL) {

    ## bring 'predictions' and 'labels' into list format,
    ## each list entry representing one x-validation run

    ## convert predictions into canonical list format
    if (is.data.frame(predictions)) {
        names(predictions) <- c()
        predictions <- as.list(predictions)
    } else if (is.matrix(predictions)) {
        predictions <- as.list(data.frame(predictions))
        names(predictions) <- c()
    } else if (is.vector(predictions) && !is.list(predictions)) {
        predictions <- list(predictions)
    } else if (!is.list(predictions)) {
        stop("Format of predictions is invalid.")
    } 
    ## if predictions is a list -> keep unaltered
  
    ## convert labels into canonical list format
    if (is.data.frame(labels)) {
        names(labels) <- c()
        labels <- as.list( labels)
    } else if (is.matrix(labels)) {
        labels <- as.list( data.frame( labels))
        names(labels) <- c()
    } else if ((is.vector(labels) ||
                is.ordered(labels) ||
                is.factor(labels)) &&
               !is.list(labels)) {
        labels <- list( labels)
    } else if (!is.list(labels)) {
        stop("Format of labels is invalid.")
    }
    ## if labels is a list -> keep unaltered

    ## Length consistency checks
    if (length(predictions) != length(labels))
      stop(paste("Number of cross-validation runs must be equal",
                 "for predictions and labels."))
    if (! all(sapply(predictions, length) == sapply(labels, length)))
      stop(paste("Number of predictions in each run must be equal",
                 "to the number of labels for each run."))
    
    ## only keep prediction/label pairs that are finite numbers
    for (i in 1:length(predictions)) {
        finite.bool <- is.finite( predictions[[i]] )
        predictions[[i]] <- predictions[[i]][ finite.bool ]
        labels[[i]] <- labels[[i]][ finite.bool ]
    }

    ## abort if 'labels' format is inconsistent across
    ## different cross-validation runs
    label.format=""  ## one of 'normal','factor','ordered'
    if (all(sapply( labels, is.factor)) &&
        !any(sapply(labels, is.ordered))) {
        label.format <- "factor"
    } else if (all(sapply( labels, is.ordered))) {
        label.format <- "ordered"
    } else if (all(sapply( labels, is.character)) || 
               all(sapply( labels, is.numeric)) ||
               all(sapply( labels, is.logical))) {
        label.format <- "normal"
    } else {
        stop(paste("Inconsistent label data type across different",
                   "cross-validation runs."))
    }
    
    ## abort if levels are not consistent across different
    ## cross-validation runs
    if (! all(sapply(labels, levels)==levels(labels[[1]])) ) {
        stop(paste("Inconsistent factor levels across different",
                   "cross-validation runs."))
    }
        
    ## convert 'labels' into ordered factors, aborting if the number
    ## of classes is not equal to 2.
    levels <- c()
    if ( label.format == "ordered" ) {
        if (!is.null(label.ordering)) {
            stop(paste("'labels' is already ordered. No additional",
                       "'label.ordering' must be supplied."))
        } else {
            levels <- levels(labels[[1]])
        }
    } else {
        if ( is.null( label.ordering )) {
            if ( label.format == "factor" ) levels <- sort(levels(labels[[1]]))
            else levels <- sort( unique( unlist( labels)))
        } else {
          ## if (!setequal( levels, label.ordering)) {
          if (!setequal( unique(unlist(labels)), label.ordering )) {
            stop("Label ordering does not match class labels.")
          }
          levels <- label.ordering
        }
        for (i in 1:length(labels)) {
            if (is.factor(labels))
              labels[[i]] <- ordered(as.character(labels[[i]]),
                                     levels=levels)
            else labels[[i]] <- ordered( labels[[i]], levels=levels)
        }

    }

    if (length(levels) != 2) {
        message <- paste("Number of classes is not equal to 2.\n",
                         "ROCR currently supports only evaluation of ",
                         "binary classification tasks.",sep="")
        stop(message)
    }

    ## determine whether predictions are continuous or categorical
    ## (in the latter case stop; scheduled for the next ROCR version)
    if (!is.numeric( unlist( predictions ))) {
        stop("Currently, only continuous predictions are supported by ROCR.")
    }

    ## compute cutoff/fp/tp data

    cutoffs <- list()
    fp <- list()
    tp <- list()
    fn <- list()
    tn <- list()
    n.pos <- list()
    n.neg <- list()
    n.pos.pred <- list()
    n.neg.pred <- list()
    for (i in 1:length(predictions)) {
        n.pos <- c( n.pos, sum( labels[[i]] == levels[2] ))
        n.neg <- c( n.neg, sum( labels[[i]] == levels[1] ))
        ans <- .compute.unnormalized.roc.curve( predictions[[i]], labels[[i]] )
        cutoffs <- c( cutoffs, list( ans$cutoffs ))
        fp <- c( fp, list( ans$fp ))
        tp <- c( tp, list( ans$tp ))
        fn <- c( fn, list( n.pos[[i]] - tp[[i]] ))
        tn <- c( tn, list( n.neg[[i]] - fp[[i]] ))
        n.pos.pred <- c(n.pos.pred, list(tp[[i]] + fp[[i]]) )
        n.neg.pred <- c(n.neg.pred, list(tn[[i]] + fn[[i]]) )
    }


    return( new("prediction", predictions=predictions,
                labels=labels,
                cutoffs=cutoffs,
                fp=fp,
                tp=tp,
                fn=fn,
                tn=tn,
                n.pos=n.pos,
                n.neg=n.neg,
                n.pos.pred=n.pos.pred,
                n.neg.pred=n.neg.pred))
}


## fast fp/tp computation based on cumulative summing
.compute.unnormalized.roc.curve <- function( predictions, labels ) {
    ## determine the labels that are used for the pos. resp. neg. class :
    pos.label <- levels(labels)[2]
    neg.label <- levels(labels)[1]

    pred.order <- order(predictions, decreasing=TRUE)
    predictions.sorted <- predictions[pred.order]
    tp <- cumsum(labels[pred.order]==pos.label)
    fp <- cumsum(labels[pred.order]==neg.label)

    ## remove fp & tp for duplicated predictions
    ## as duplicated keeps the first occurrence, but we want the last, two
    ## rev are used.
    ## Highest cutoff (Infinity) corresponds to tp=0, fp=0
    dups <- rev(duplicated(rev(predictions.sorted)))
    tp <- c(0, tp[!dups])
    fp <- c(0, fp[!dups])
    cutoffs <- c(Inf, predictions.sorted[!dups])
    
    return(list( cutoffs=cutoffs, fp=fp, tp=tp ))
}

## ------------------------------------------------------------------------
## classical machine learning contingency table measures
## ------------------------------------------------------------------------

.performance.accuracy <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    
      list( cutoffs, (tn+tp) / length(predictions) )
  }

.performance.error.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

      list( cutoffs, (fn+fp) / length(predictions) )
  }

.performance.false.positive.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    
      list( cutoffs, fp / n.neg )
  }

.performance.true.positive.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    
      list( cutoffs, tp / n.pos )
  }

.performance.false.negative.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    
      list( cutoffs, fn / n.pos )
  }

.performance.true.negative.rate <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    
      list( cutoffs, tn / n.neg )
  }

.performance.positive.predictive.value <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    
      ppv <- tp / (fp + tp)
      list( cutoffs, ppv )
  }

.performance.negative.predictive.value <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
    
      npv <- tn / (tn + fn)
      list( cutoffs, npv )
  }

.performance.prediction.conditioned.fallout <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
      ppv <- .performance.positive.predictive.value(predictions, labels,
                                                    cutoffs, fp, tp, fn, tn,
                                                    n.pos, n.neg, n.pos.pred,
                                                    n.neg.pred)[[2]]
      list( cutoffs, 1 - ppv )
  }

.performance.prediction.conditioned.miss <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
      npv <- .performance.negative.predictive.value(predictions, labels,
                                                    cutoffs, fp, tp, fn, tn,
                                                    n.pos, n.neg, n.pos.pred,
                                                    n.neg.pred)[[2]]
      list( cutoffs, 1 - npv )
  }

## ------------------------------------------------------------------------
## ...not actually performance measures, but very useful as a second axis
## against which to plot a "real" performance measure
## (popular example: lift charts)
## ------------------------------------------------------------------------

.performance.rate.of.positive.predictions <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

      list( cutoffs, n.pos.pred / (n.pos + n.neg) )
  }

.performance.rate.of.negative.predictions <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

      list( cutoffs, n.neg.pred / (n.pos + n.neg) )
  }


## ------------------------------------------------------------------------
## Classical statistical contingency table measures
## ------------------------------------------------------------------------

.performance.phi <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

      list(cutoffs,
           (tn*tp - fn*fp) / sqrt(n.pos * n.neg * n.pos.pred * n.neg.pred) )
  }

.performance.mutual.information <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

      n.samples <- n.pos + n.neg
      mi <- c()
      for (k in 1:length(cutoffs)) {
          kij <- rbind( c(tn[k],fn[k]), c(fp[k],tp[k]) )

          ki.j. <- rbind(c(n.neg * n.neg.pred[k], n.neg.pred[k] * n.pos),
                         c(n.neg * n.pos.pred[k], n.pos * n.pos.pred[k]))

          log.matrix <- log2( kij / ki.j.)
          log.matrix[kij/ki.j.==0] <- 0
          
          mi <- c(mi,  log2(n.samples) + sum( kij * log.matrix) / n.samples  )
      }

      list( cutoffs, mi )
  }


.performance.chisq <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

      chisq <- c()
      for (i in 1:length(cutoffs)) {
          A <- rbind( c( tn[i], fn[i]), c(fp[i], tp[i]) )
          chisq <- c(chisq, chisq.test(A, correct=FALSE)$statistic )
      }
      list( cutoffs, chisq )
  }

.performance.odds.ratio <- 
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
      

    list( cutoffs, tp * tn / (fn * fp) )
    
}

## ------------------------------------------------------------------------
## Other measures based on contingency tables
## ------------------------------------------------------------------------

.performance.lift <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
      
      n.samples <- n.pos + n.neg
      list( cutoffs, (tp / n.pos) / (n.pos.pred / n.samples) )
  }

.performance.f <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, alpha) {

      prec <- .performance.positive.predictive.value(predictions, labels,
                                                     cutoffs, fp, tp, fn, tn,
                                                     n.pos, n.neg, n.pos.pred,
                                                     n.neg.pred)[[2]]
      list( cutoffs,  1/ ( alpha*(1/prec) + (1-alpha)*(1/(tp/n.pos))  ) )
  }

.performance.rocconvexhull <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
      
      x <- fp / n.neg
      y <- tp / n.pos

      finite.bool <- is.finite(x) & is.finite(y)
      x <- x[ finite.bool ]
      y <- y[ finite.bool ]
      if (length(x) < 2) {
          stop("Not enough distinct predictions to compute ROC convex hull.")
      }

      ## keep only points on the convex hull
      ind <- chull(x, y)
      x.ch <- x[ind]
      y.ch <- y[ind]

      ## keep only convex hull points above the diagonal, except (0,0)
      ## and (1,1)
      ind.upper.triangle <- x.ch < y.ch
      x.ch <- c(0, x.ch[ind.upper.triangle], 1)
      y.ch <- c(0, y.ch[ind.upper.triangle], 1)

      ## sort remaining points by ascending x value
      ind <- order(x.ch)
      x.ch <- x.ch[ind]
      y.ch <- y.ch[ind]

      list( x.ch, y.ch )
  }

## ----------------------------------------------------------------------------
## Cutoff-independent measures
## ----------------------------------------------------------------------------

.performance.auc <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, fpr.stop) {
      
      x <- fp / n.neg
      y <- tp / n.pos

      finite.bool <- is.finite(x) & is.finite(y)
      x <- x[ finite.bool ]
      y <- y[ finite.bool ]
      if (length(x) < 2) {
          stop(paste("Not enough distinct predictions to compute area",
                     "under the ROC curve."))
      }

      if (fpr.stop < 1) {
        ind <- max(which( x <= fpr.stop ))
        tpr.stop <- approxfun( x[ind:(ind+1)], y[ind:(ind+1)] )(fpr.stop)
        x <- c(x[1:ind], fpr.stop)
        y <- c(y[1:ind], tpr.stop)
      }
      
      ans <- list()
      auc <- 0
      for (i in 2:length(x)) {
          auc <- auc + 0.5 * (x[i] - x[i-1]) * (y[i] + y[i-1])
      }
      ans <- list( c(), auc)
      names(ans) <- c("x.values","y.values")
      return(ans)
  }

.performance.precision.recall.break.even.point <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

      pred <- prediction( predictions, labels)
      perf <- performance( pred, measure="prec", x.measure="rec")
      x <- rev(perf@x.values[[1]])
      y <- rev(perf@y.values[[1]])
      alpha <- rev(perf@alpha.values[[1]])

      finite.bool <- is.finite(alpha) & is.finite(x) & is.finite(y)
      x <- x[ finite.bool ]
      y <- y[ finite.bool ]
      alpha <- alpha[ finite.bool ]

      if (length(x) < 2) {
          stop(paste("Not enough distinct predictions to compute",
                     "precision/recall intersections."))
      }
      intersection.cutoff <- c()
      intersection.pr <- c()
      
      ## find all intersection points by looking at all intervals (i,i+1):
      ## if the difference function between x and y has different signs at the
      ## interval boundaries, then an intersection point is in the interval;
      ## compute as the root of the difference function
      if ( (x[1]-y[1]) == 0) {
          intersection.cutoff <- c( alpha[1] )
          intersection.pr <- c( x[1] )
      }

      for (i in (1:(length(alpha)-1))) {
          if ((x[i+1]-y[i+1]) == 0) {
              intersection.cutoff <- c( intersection.cutoff, alpha[i+1] )
              intersection.pr <- c( intersection.pr, x[i+1] )
          } else if ((x[i]-y[i])*(x[i+1]-y[i+1]) < 0 ) {
              ans <- uniroot(approxfun(c(alpha[i], alpha[i+1] ),
                                       c(x[i]-y[i], x[i+1]-y[i+1])),
                             c(alpha[i],alpha[i+1]))
              intersection.cutoff <- c(intersection.cutoff, ans$root)
              intersection.pr <- c(intersection.pr, ans$f.root)
          }
      }

      list( rev(intersection.cutoff), rev(intersection.pr) )
  }


.performance.calibration.error <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, window.size) {

      if (window.size > length(predictions)) {
          stop("Window size exceeds number of predictions.")
      }
      if (min(predictions)<0 || max(predictions)>1) {
          stop("Calibration error needs predictions between 0 and 1")
      }
      
      pos.label <- levels(labels)[2]
      neg.label <- levels(labels)[1]

      ordering <- rev(order( predictions ))
      predictions <- predictions[ ordering ]
      labels <- labels[ ordering ]

      median.cutoffs <- c()
      calibration.errors <- c()

      for (left.index in 1 : (length(predictions) - window.size+1) ) {
          right.index <- left.index + window.size - 1
          pos.fraction <-
            sum(labels[left.index : right.index] == pos.label) / window.size
          mean.prediction <- mean( predictions[ left.index : right.index ] )

          calibration.errors <- c(calibration.errors,
                                  abs(pos.fraction - mean.prediction))
          median.cutoffs <- c(median.cutoffs,
                              median(predictions[left.index:right.index]))
      }
      list( median.cutoffs, calibration.errors )
  }


.performance.mean.cross.entropy <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

      if (! all(levels(labels)==c(0,1)) ||
          any(predictions<0) || any(predictions>1) ) {
          stop(paste("Class labels need to be 0 and 1 and predictions between",
                     "0 and 1 for mean cross entropy."))
      }
      
      pos.label <- levels(labels)[2]
      neg.label <- levels(labels)[1]
    
      list( c(), - 1/length(predictions) *
           (sum( log( predictions[which(labels==pos.label)] ))  +
            sum( log( 1 - predictions[which(labels==neg.label)] ))) )
  }


.performance.root.mean.squared.error <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
      ## convert labels from factor to numeric values
      labels <- as.numeric(levels(labels))[labels]
      if (any(is.na(labels))) {
          stop("For rmse predictions have to be numeric.")
      }
      list( c(),  sqrt( 1/length(predictions) *
                      sum( (predictions - labels)^2 ))  )
  }

## ----------------------------------------------------------------------------
## Derived measures:
## ----------------------------------------------------------------------------

.performance.sar <- function( predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {

    pred <- prediction( predictions, labels)
    perf.acc <- performance( pred, measure="acc")
    perf.rmse <- performance( pred, measure="rmse")
    perf.auc <- performance( pred, measure="auc")

    list(cutoffs,
         1/3 * (perf.acc@y.values[[1]] +
                (1 - perf.rmse@y.values[[1]]) +
                perf.auc@y.values[[1]]))
}

## ----------------------------------------------------------------------------
## Measures taking into account actual cost considerations
## ----------------------------------------------------------------------------

.performance.expected.cost <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred) {
      
      ## kick out suboptimal values (i.e. fpr/tpr pair for which another one
      ## with same fpr and higher tpr exists, 
      ## or one for which one with same tpr but lower fpr exists

      if (n.neg==0 || n.pos==0) {
          stop(paste("At least one positive and one negative sample are",
                     "needed to compute a cost curve."))
      }
      fpr <- fp / n.neg
      tpr <- tp / n.pos

      ## sort by fpr (ascending), in case of ties by descending tpr
      ind <- order(fpr,-tpr)
      
      fpr <- fpr[ind]
      tpr <- tpr[ind]
      ## for tied fprs, only the one with the highest tpr is kept
      ind <- !duplicated(fpr)
      fpr <- fpr[ind]
      tpr <- tpr[ind]

      ## for tied tprs, only keep the one with the lowest fpr
      ind <- order(-tpr,fpr)
      fpr <- fpr[ind]
      tpr <- tpr[ind]
      ind <- !duplicated(tpr)
      fpr <- fpr[ind]
      tpr <- tpr[ind]

      if (!any(0==fpr & 0==tpr)) {
          fpr <- c(0,fpr)
          tpr <- c(0,tpr)
      }
      if (!any(1==fpr & 1==tpr)) {
          fpr <- c(fpr,1)
          tpr <- c(tpr,1)
      }
      
      ## compute all functions
      f <- list()
      for (i in 1:length(fpr)) {
          f <- c(f, .construct.linefunct( 0, fpr[i], 1, 1-tpr[i] ))
      }
      
      ## compute all intersection points
      x.values <- c()
      y.values <- c()
      for (i in 1:(length(fpr)-1)) {
          for (j in (i+1):length(fpr)) {
              ans <- .intersection.point( f[[i]], f[[j]] )
              if (all(is.finite(ans))) {
                  y.values.at.current.x <- c()
                  for (k in 1:length(f)) {
                      y.values.at.current.x <- c(y.values.at.current.x,
                                                 f[[k]](ans[1]))
                  }
                  if (abs(ans[2] - min(y.values.at.current.x )) <
                      sqrt(.Machine$double.eps)) {
                      
                      x.values <- c(x.values, ans[1])
                      y.values <- c(y.values, ans[2])
                  }
              }
          }
      }

      if (!any(0==x.values & 0==y.values)) {
          x.values <- c(0,x.values)
          y.values <- c(0,y.values)
      }
      if (!any(1==x.values & 0==y.values)) {
          x.values <- c(x.values,1)
          y.values <- c(y.values,0)
      }

      ind <- order( x.values)
      list( x.values[ind], y.values[ind] )
  }


.performance.cost <-
  function(predictions, labels, cutoffs, fp, tp, fn, tn,
           n.pos, n.neg, n.pos.pred, n.neg.pred, cost.fp, cost.fn) {
      
    n.samples <- n.pos + n.neg
    cost <- ((n.pos / n.samples) * (fn / n.pos) * cost.fn +
             (n.neg / n.samples) * (fp / n.neg) * cost.fp)
    list( cutoffs, cost )
}




setClass("prediction",
         representation(predictions = "list",
                        labels      = "list",
                        cutoffs     = "list",
                        fp          = "list",
                        tp          = "list",
                        tn          = "list",
                        fn          = "list",
                        n.pos       = "list",
                        n.neg       = "list",
                        n.pos.pred  = "list",
                        n.neg.pred  = "list"))

setClass("performance",
         representation(x.name       = "character",
                        y.name       = "character",
                        alpha.name   = "character",
                        x.values     = "list",
                        y.values     = "list",
                        alpha.values = "list" ))



## ---------------------------------------------------------------------------
## Dealing with argument lists, especially '...'
## ---------------------------------------------------------------------------

## return list of selected arguments, skipping those that
## are not present in arglist
.select.args <- function( arglist, args.to.select, complement=FALSE) {
    match.bool <- names(arglist) %in% args.to.select
    if (complement==TRUE) match.bool <- !match.bool
    return( arglist[ match.bool] )
}

## return arguments in arglist which match prefix, with prefix removed
## ASSUMPTION: prefix is separated from rest by a '.'; this is removed along
## with the prefix
.select.prefix <- function( arglist, prefixes, complement=FALSE ) {
    match.expr <- paste(paste('(^',prefixes,'\\.)',sep=""),collapse='|')
    match.bool <- (1:length(arglist)) %in% grep( match.expr, names(arglist) )
    if (complement==TRUE) match.bool <- !match.bool
    arglist <- arglist[ match.bool]
    names(arglist) <- sub( match.expr, '', names(arglist))
    
    return( arglist )
}

.garg <- function( arglist, arg, i=1) {
    if (is.list(arglist[[arg]])) arglist[[ arg ]][[i]]
    else arglist[[ arg ]]
}

.sarg <- function( arglist, ...) {
    ll <- list(...)
    for (argname in names(ll) ) {
        arglist[[ argname ]] <- ll[[ argname ]]
    }
    return(arglist)
}

.farg <- function( arglist, ...) {
    ll <- list(...)
    for (argname in names(ll) ) {
        if (length(arglist[[argname]])==0)
          arglist[[ argname ]] <- ll[[ argname ]]
    }
    return(arglist)
}

.slice.run <- function( arglist, runi=1) {
    r <- lapply( names(arglist), function(name) .garg( arglist, name, runi))
    names(r) <- names(arglist)
    r
}

performance <- function(prediction.obj, measure,
                        x.measure="cutoff", ...) {

    ## define the needed environments
    envir.list <- .define.environments()
    long.unit.names <- envir.list$long.unit.names
    function.names <- envir.list$function.names
    obligatory.x.axis <- envir.list$obligatory.x.axis
    optional.arguments <- envir.list$optional.arguments
    default.values <- envir.list$default.values
    
    
    
    ## abort, if attempt is made to use a measure that has an obligatory
    ## x.axis as the x.measure (cannot be combined)
    if (exists( x.measure, where=obligatory.x.axis, inherits=FALSE )) {
        message <- paste("The performance measure",
                         x.measure,
                         "can only be used as 'measure', because it has",
                         "the following obligatory 'x.measure':\n",
                         get( x.measure, envir=obligatory.x.axis))
        stop(message)
    }

    ## if measure is a performance measure with obligatory x.axis, then
    ## enforce this axis:
    if (exists( measure, where=obligatory.x.axis, inherits=FALSE )) {
        x.measure <- get( measure, envir=obligatory.x.axis )
    }

    if (x.measure == "cutoff" ||
        exists( measure, where=obligatory.x.axis, inherits=FALSE )) {

        ## fetch from '...' any optional arguments for the performance
        ## measure at hand that are given, otherwise fill up the default values
        optional.args <- list(...)
        argnames <- c()
        if ( exists( measure, where=optional.arguments, inherits=FALSE )) {
            argnames <- get( measure, envir=optional.arguments )
            default.arglist <- list()
            for (i in 1:length(argnames)) {
                default.arglist <- c(default.arglist,
                                     get(paste(measure,":",argnames[i],sep=""),
                                         envir=default.values, inherits=FALSE))
            }
            names(default.arglist) <- argnames

            for (i in 1:length(argnames)) {
                templist <- list(optional.args,
                                 default.arglist[[i]])
                names(templist) <- c('arglist', argnames[i])
                
                optional.args <- do.call('.farg', templist)
            }
        }
        optional.args <- .select.args( optional.args, argnames )
        
        ## determine function name
        function.name <- get( measure, envir=function.names )

        ## for each x-validation run, compute the requested performance measure
        x.values <- list()
        y.values <- list()
        for (i in 1:length( prediction.obj@predictions )) {
            argumentlist <- .sarg(optional.args,
                                  predictions= prediction.obj@predictions[[i]],
                                  labels= prediction.obj@labels[[i]],
                                  cutoffs= prediction.obj@cutoffs[[i]],
                                  fp= prediction.obj@fp[[i]],
                                  tp= prediction.obj@tp[[i]],
                                  fn= prediction.obj@fn[[i]],
                                  tn= prediction.obj@tn[[i]],
                                  n.pos= prediction.obj@n.pos[[i]],
                                  n.neg= prediction.obj@n.neg[[i]],
                                  n.pos.pred= prediction.obj@n.pos.pred[[i]],
                                  n.neg.pred= prediction.obj@n.neg.pred[[i]])

            ans <- do.call( function.name, argumentlist )

            if (!is.null(ans[[1]])) x.values <- c( x.values, list( ans[[1]] ))
            y.values <- c( y.values, list( ans[[2]] ))
        }

        if (! (length(x.values)==0 || length(x.values)==length(y.values)) ) {
            stop("Consistency error.")
        }
        
        ## create a new performance object
        return( new("performance",
                    x.name       = get( x.measure, envir=long.unit.names ),
                    y.name       = get( measure, envir=long.unit.names ),
                    alpha.name   = "none",
                    x.values     = x.values,
                    y.values     = y.values,
                    alpha.values = list() ))
    } else {
        perf.obj.1 <- performance( prediction.obj, measure=x.measure, ... )
        perf.obj.2 <- performance( prediction.obj, measure=measure, ... )
        return( .combine.performance.objects( perf.obj.1, perf.obj.2 ) )
    }
}

.combine.performance.objects <- function( p.obj.1, p.obj.2 ) {
    ## some checks for misusage (in any way, this function is
    ## only for internal use)
    if ( p.obj.1@x.name != p.obj.2@x.name ) {
        stop("Error: Objects need to have identical x axis.")
    }
    if ( p.obj.1@alpha.name != "none" || p.obj.2@alpha.name != "none") {
        stop("Error: At least one of the two objects has already been merged.")
    }
    if (length(p.obj.1@x.values) != length(p.obj.2@x.values)) {
        stop(paste("Only performance objects with identical number of",
                   "cross-validation runs can be combined."))
    }

    x.values <- list()
    x.name <- p.obj.1@y.name
    y.values <- list()
    y.name <- p.obj.2@y.name
    alpha.values <- list()
    alpha.name <- p.obj.1@x.name

    for (i in 1:length( p.obj.1@x.values )) {
        x.values.1 <- p.obj.1@x.values[[i]]
        y.values.1 <- p.obj.1@y.values[[i]]
        x.values.2 <- p.obj.2@x.values[[i]]
        y.values.2 <- p.obj.2@y.values[[i]]

        ## cutoffs of combined object = merged cutoffs of simple objects
        cutoffs <- sort( unique( c(x.values.1, x.values.2)), decreasing=TRUE )

        ## calculate y.values at cutoffs using step function
        y.values.int.1 <- approxfun(x.values.1, y.values.1,
                                    method="constant",f=1,rule=2)(cutoffs)
        y.values.int.2 <- approxfun(x.values.2, y.values.2,
                                    method="constant",f=1,rule=2)(cutoffs)

        ## 'approxfun' ignores NA and NaN
        objs <- list( y.values.int.1, y.values.int.2)
        objs.x <- list( x.values.1, x.values.2 )
        na.cutoffs.1.bool <- is.na( y.values.1) & !is.nan( y.values.1 )
        nan.cutoffs.1.bool <- is.nan( y.values.1)
        na.cutoffs.2.bool <- is.na( y.values.2) & !is.nan( y.values.2 )
        nan.cutoffs.2.bool <- is.nan( y.values.2)
        bools <- list(na.cutoffs.1.bool, nan.cutoffs.1.bool,
                      na.cutoffs.2.bool, nan.cutoffs.2.bool)
        values <- c(NA,NaN,NA,NaN)
        
        for (j in 1:4) {
            for (k in which(bools[[j]])) {
                interval.max <- objs.x[[ ceiling(j/2) ]][k]
                interval.min <- -Inf
                if (k < length(objs.x[[ ceiling(j/2) ]])) {
                    interval.min <- objs.x[[ ceiling(j/2) ]][k+1]
                }
                objs[[ ceiling(j/2) ]][cutoffs <= interval.max &
                                       cutoffs > interval.min ] <- values[j]
            }
        }

        alpha.values <- c(alpha.values, list(cutoffs))
        x.values <- c(x.values, list(objs[[1]]))
        y.values <- c(y.values, list(objs[[2]]))
    }
    
    return( new("performance",
                x.name=x.name, y.name=y.name,
                alpha.name=alpha.name, x.values=x.values,
                y.values=y.values, alpha.values=alpha.values))
}

.define.environments <- function() {
    ## There are five environments: long.unit.names, function.names,
    ## obligatory.x.axis, optional.arguments, default.values
    
    ## Define long names corresponding to the measure abbreviations.
    long.unit.names <- new.env()
    assign("none","None", envir=long.unit.names)
    assign("cutoff", "Cutoff", envir=long.unit.names)
    assign("acc", "Accuracy", envir=long.unit.names)
    assign("err", "Error Rate", envir=long.unit.names)
    assign("fpr", "False positive rate", envir=long.unit.names)
    assign("tpr", "True positive rate", envir=long.unit.names)
    assign("rec", "Recall", envir=long.unit.names)
    assign("sens", "Sensitivity", envir=long.unit.names)
    assign("fnr", "False negative rate", envir=long.unit.names)
    assign("tnr", "True negative rate", envir=long.unit.names)
    assign("spec", "Specificity", envir=long.unit.names)
    assign("ppv", "Positive predictive value", envir=long.unit.names)
    assign("prec", "Precision", envir=long.unit.names)
    assign("npv", "Negative predictive value", envir=long.unit.names)
    assign("fall", "Fallout", envir=long.unit.names)
    assign("miss", "Miss", envir=long.unit.names)
    assign("pcfall", "Prediction-conditioned fallout", envir=long.unit.names)
    assign("pcmiss", "Prediction-conditioned miss", envir=long.unit.names)
    assign("rpp", "Rate of positive predictions", envir=long.unit.names)
    assign("rnp", "Rate of negative predictions", envir=long.unit.names)
    assign("auc","Area under the ROC curve", envir=long.unit.names)
    assign("cal", "Calibration error", envir=long.unit.names)
    assign("mwp", "Median window position", envir=long.unit.names)
    assign("prbe","Precision/recall break-even point", envir=long.unit.names)
    assign("rch", "ROC convex hull", envir=long.unit.names)
    assign("mxe", "Mean cross-entropy", envir=long.unit.names)
    assign("rmse","Root-mean-square error", envir=long.unit.names)
    assign("phi", "Phi correlation coefficient", envir=long.unit.names)
    assign("mat","Matthews correlation coefficient", envir=long.unit.names)
    assign("mi", "Mutual information", envir=long.unit.names)
    assign("chisq", "Chi-square test statistic", envir=long.unit.names)
    assign("odds","Odds ratio", envir=long.unit.names)
    assign("lift", "Lift value", envir=long.unit.names)
    assign("f","Precision-Recall F measure", envir=long.unit.names)
    assign("sar", "SAR", envir=long.unit.names)
    assign("ecost", "Expected cost", envir=long.unit.names)
    assign("cost", "Explicit cost", envir=long.unit.names)

    ## Define function names corresponding to the measure abbreviations.
    function.names <- new.env()
    assign("acc", ".performance.accuracy", envir=function.names)
    assign("err", ".performance.error.rate", envir=function.names)
    assign("fpr", ".performance.false.positive.rate", envir=function.names)
    assign("tpr", ".performance.true.positive.rate", envir=function.names)
    assign("rec", ".performance.true.positive.rate", envir=function.names)
    assign("sens", ".performance.true.positive.rate", envir=function.names)
    assign("fnr", ".performance.false.negative.rate", envir=function.names)
    assign("tnr", ".performance.true.negative.rate", envir=function.names)
    assign("spec", ".performance.true.negative.rate", envir=function.names)
    assign("ppv", ".performance.positive.predictive.value",
           envir=function.names)
    assign("prec", ".performance.positive.predictive.value",
           envir=function.names)
    assign("npv", ".performance.negative.predictive.value",
           envir=function.names)
    assign("fall", ".performance.false.positive.rate", envir=function.names)
    assign("miss", ".performance.false.negative.rate", envir=function.names)
    assign("pcfall", ".performance.prediction.conditioned.fallout",
           envir=function.names)
    assign("pcmiss", ".performance.prediction.conditioned.miss",
           envir=function.names)
    assign("rpp", ".performance.rate.of.positive.predictions",
           envir=function.names)
    assign("rnp", ".performance.rate.of.negative.predictions",
           envir=function.names)
    assign("auc", ".performance.auc", envir=function.names)
    assign("cal", ".performance.calibration.error", envir=function.names)
    assign("prbe", ".performance.precision.recall.break.even.point",
           envir=function.names)
    assign("rch", ".performance.rocconvexhull", envir=function.names)
    assign("mxe", ".performance.mean.cross.entropy", envir=function.names)
    assign("rmse", ".performance.root.mean.squared.error",
           envir=function.names)
    assign("phi", ".performance.phi", envir=function.names)
    assign("mat", ".performance.phi", envir=function.names)
    assign("mi", ".performance.mutual.information", envir=function.names)
    assign("chisq", ".performance.chisq", envir=function.names)
    assign("odds", ".performance.odds.ratio", envir=function.names)
    assign("lift", ".performance.lift", envir=function.names)
    assign("f", ".performance.f", envir=function.names)
    assign("sar", ".performance.sar", envir=function.names)
    assign("ecost", ".performance.expected.cost", envir=function.names)
    assign("cost", ".performance.cost", envir=function.names)

    ## If a measure comes along with an obligatory x axis (including "none"),
    ## list it here.
    obligatory.x.axis <- new.env()
    assign("mxe", "none", envir=obligatory.x.axis)
    assign("rmse", "none", envir=obligatory.x.axis)
    assign("prbe", "none", envir=obligatory.x.axis)
    assign("auc", "none", envir=obligatory.x.axis)
    assign("rch","none", envir=obligatory.x.axis)
    ## ecost requires probability cost function as x axis, which is handled
    ## implicitly, not as an explicit performance measure.
    assign("ecost","none", envir=obligatory.x.axis)  
    
    ## If a measure has optional arguments, list the names of the
    ## arguments here.
    optional.arguments <- new.env()
    assign("cal", "window.size", envir=optional.arguments)
    assign("f", "alpha", envir=optional.arguments)
    assign("cost", c("cost.fp", "cost.fn"), envir=optional.arguments)
    assign("auc", "fpr.stop", envir=optional.arguments)
        
    ## If a measure has additional arguments, list the default values
    ## for them here. Naming convention: e.g. "cal" has an optional
    ## argument "window.size" the key to use here is "cal:window.size"
    ## (colon as separator)
    default.values <- new.env()
    assign("cal:window.size", 100, envir=default.values)
    assign("f:alpha", 0.5, envir=default.values)
    assign("cost:cost.fp", 1, envir=default.values)
    assign("cost:cost.fn", 1, envir=default.values)
    assign("auc:fpr.stop", 1, envir=default.values) 
    
    list(long.unit.names=long.unit.names, function.names=function.names,
         obligatory.x.axis=obligatory.x.axis,
         optional.arguments=optional.arguments,
         default.values=default.values)
}
```

```{r ex="create_d", type="sample-code"}
# df, including pred, is loaded

# Use the prediction function from RCOR (preloaded) to make a prediction object
# The first argument is our prediction from the regression, and the second is
# The actual outcome
ROCpred <- prediction(, )

# Use the performance function (preloaded) from ROCR along with ROCpred
# The created object will be used to plot a ROC curve, and thus needs
# true positive rate ('tpr') and false positive rate ('fpr')
ROCperf <- performance(, '','')

# Create a data frame for passing to ggplot2 for plotting the curve
# This part is already done
df_ROC <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                     TruePositive=c(ROCperf@y.values[[1]]))

# Calculate ROC AUC ('auc') using performance() and ROCpred
AUC <- performance(, '')

# Display the AUC value -- this line is already done
# The @ is in place of a $, due to ROCR using S4 objects instead of lists
AUC@y.values[[1]]
#END
```

```{r ex="create_d", type="solution"}
# df, including pred, is loaded

# Use the prediction function from RCOR (preloaded) to make a prediction object
# The first argument is our prediction from the regression, and the second is
# The actual outcome
ROCpred <- prediction(as.numeric(df$pred), as.numeric(df$revtq_down))

# Use the performance function (preloaded) from ROCR along with ROCpred
# The created object will be used to plot a ROC curve, and thus needs
# true positive rate ('tpr') and false positive rate ('fpr')
ROCperf <- performance(ROCpred, 'tpr','fpr')

# Create a data frame for passing to ggplot2 for plotting the curve
# This part is already done
df_ROC <- data.frame(FalsePositive=c(ROCperf@x.values[[1]]),
                     TruePositive=c(ROCperf@y.values[[1]]))

# Calculate ROC AUC ('auc') using performance() and ROCpred
AUC <- performance(ROCpred, 'auc')

# Display the AUC value -- this line is already done
# The @ is in place of a $, due to ROCR using S4 objects instead of lists
AUC@y.values[[1]]
#END
```

```{r ex="create_d", type="sct"}
# Template based on https://www.rdocumentation.org/packages/testwhat/versions/4.1.1
# Check if something is explicitly typed

test_expression_output("ROCpred", incorrect_msg="Your regression model isn't quite right.")
test_expression_output("ROCperf", incorrect_msg="Your regression model isn't quite right.")
test_expression_output("AUC@y.values[[1]]", incorrect_msg="Your regression model isn't quite right.")


#test_student_typed('revtq_down ~ recession_lag',not_typed_msg='Make sure to regress `revtq_down` on `recession_lag`.')


# test_student_typed('x <- 2', not_typed_msg='')

# Check if function was used in input code
# test_function('c',incorrect_msg='')  

# Requires an object `x` to have the same value as the solution
# test_object("x",incorrect_msg = "",undefined_msg = "")  

# Requires an onject with the same value of `x` in the solution
# test_an_object("x",undefined_msg="")

# Checks if output of student's code contains given evaluated expression
# test_output_contains("x",incorrect_msg = "")

# Check if a vector of predefined objects are unchanged
# test_predefined_objects(c('x','y'),incorrect_msg="Don't onverwrite the predefined variables")

# Checks for a regex pattern in trhe output
# test_output_regex(pattern,fixed=F, times=1, incorrect_msg='')

# Can check an arbitrary expression across both solution and student code
#test_expression_output("typeof(company_name)", incorrect_msg="Did you store textual data in `company_name`?")

test_error()
success_msg("Awesome!")

# Other functions to note:
#     - test_or(a,b) -- checks if either test a or test b pass
#     - test_ggplot() -- can check if plots are correct
#     - test_function() -- can also check included parameters
#     - test_loop() -- checking for and while loops
#     - test_library_function('package', not_called_msg='',incorrect_msg='')
#     - test_if_else() -- checking if statements
#     - test_expression_error() -- can check if functions are properly defined
#     - test_operator('operator',), not_called_msg='',incorrect_msg='')
#     - test_function_definition() -- rigorously check defined function
#     - test_data_frame() -- check if dataframe [columns] are equivalent
#     - test_function_result, test_expression_result
```

### Exercise 4: Plotting a ROC curve

```{r ex="create_e", type="pre-exercise-code"}
library(ggplot2)
df_ROC <- data.frame(FalsePositive=c(0, 0.00543478260869565, 0.0108695652173913, 0.0108695652173913, 0.0108695652173913, 0.016304347826087, 0.0217391304347826, 0.0217391304347826, 0.0271739130434783, 0.0326086956521739, 0.0380434782608696, 0.0434782608695652, 0.0489130434782609, 0.0543478260869565, 0.0543478260869565, 0.0597826086956522, 0.0597826086956522, 0.0652173913043478, 0.0706521739130435, 0.0760869565217391, 0.0815217391304348, 0.0869565217391304, 0.0923913043478261, 0.0978260869565217, 0.103260869565217, 0.108695652173913, 0.114130434782609, 0.119565217391304, 0.125, 0.130434782608696, 0.135869565217391, 0.141304347826087, 0.146739130434783, 0.152173913043478, 0.157608695652174, 0.16304347826087, 0.168478260869565, 0.173913043478261, 0.179347826086957, 0.184782608695652, 0.190217391304348, 0.195652173913043, 0.201086956521739, 0.206521739130435, 0.21195652173913, 0.217391304347826, 0.222826086956522, 0.228260869565217, 0.233695652173913, 0.239130434782609, 0.244565217391304, 0.25, 0.255434782608696, 0.260869565217391, 0.266304347826087, 0.271739130434783, 0.277173913043478, 0.282608695652174, 0.28804347826087, 0.293478260869565, 0.298913043478261, 0.304347826086957, 0.309782608695652, 0.315217391304348, 0.320652173913043, 0.326086956521739, 0.331521739130435, 0.33695652173913, 0.342391304347826, 0.347826086956522, 0.353260869565217, 0.358695652173913, 0.364130434782609, 0.369565217391304, 0.375, 0.380434782608696, 0.385869565217391, 0.391304347826087, 0.396739130434783, 0.402173913043478, 0.407608695652174, 0.41304347826087, 0.418478260869565, 0.423913043478261, 0.429347826086957, 0.434782608695652, 0.440217391304348, 0.445652173913043, 0.451086956521739, 0.456521739130435, 0.46195652173913, 0.467391304347826, 0.472826086956522, 0.478260869565217, 0.483695652173913, 0.489130434782609, 0.494565217391304, 0.5, 0.505434782608696, 0.510869565217391, 0.516304347826087, 0.521739130434783, 0.527173913043478, 0.532608695652174, 0.53804347826087, 0.543478260869565, 0.548913043478261, 0.554347826086957, 0.559782608695652, 0.565217391304348, 0.570652173913043, 0.576086956521739, 0.581521739130435, 0.58695652173913, 0.592391304347826, 0.597826086956522, 0.603260869565217, 0.608695652173913, 0.614130434782609, 0.619565217391304, 0.625, 0.630434782608696, 0.635869565217391, 0.641304347826087, 0.646739130434783, 0.652173913043478, 0.657608695652174, 0.66304347826087, 0.668478260869565, 0.673913043478261, 0.679347826086957, 0.684782608695652, 0.690217391304348, 0.695652173913043, 0.701086956521739, 0.706521739130435, 0.71195652173913, 0.717391304347826, 0.722826086956522, 0.728260869565217, 0.733695652173913, 0.739130434782609, 0.744565217391304, 0.75, 0.755434782608696, 0.760869565217391, 0.766304347826087, 0.771739130434783, 0.777173913043478, 1),
                     TruePositive=c(0, 0, 0, 0.2, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.8, 0.8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
```

```{r ex="create_e", type="sample-code"}
# df_ROC, from the previous exercise, is loaded
df_ROC[1:5,]

# Plot the ROC curve -- fill in the below command
ggplot(data=, aes(x=, y=)) +
  geom_line() +
  geom_abline(slope=1)

#END
```

```{r ex="create_e", type="solution"}
# df_ROC, from the previous exercise, is loaded
df_ROC[1:5,]

# Plot the ROC curve -- fill in the below command
ggplot(data=df_ROC, aes(x=FalsePositive, y=TruePositive)) +
  geom_line() +
  geom_abline(slope=1)

#END
```

```{r ex="create_e", type="sct"}
# Template based on https://www.rdocumentation.org/packages/testwhat/versions/4.1.1
# Check if something is explicitly typed

test_student_typed('data=df_ROC',not_typed_msg='Make sure to use `df_ROC` as the data frame.')
test_student_typed('x=FalsePositive',not_typed_msg='Make sure to put `FalsePositive` as the data plotted on the x-axis')
test_student_typed('y=FalsePositive',not_typed_msg='Make sure to put `TruePositive` as the data plotted on the y-axis')

# test_student_typed('x <- 2', not_typed_msg='')

# Check if function was used in input code
# test_function('c',incorrect_msg='')  

# Requires an object `x` to have the same value as the solution
# test_object("x",incorrect_msg = "",undefined_msg = "")  

# Requires an onject with the same value of `x` in the solution
# test_an_object("x",undefined_msg="")

# Checks if output of student's code contains given evaluated expression
# test_output_contains("x",incorrect_msg = "")

# Check if a vector of predefined objects are unchanged
# test_predefined_objects(c('x','y'),incorrect_msg="Don't onverwrite the predefined variables")

# Checks for a regex pattern in trhe output
# test_output_regex(pattern,fixed=F, times=1, incorrect_msg='')

# Can check an arbitrary expression across both solution and student code
#test_expression_output("typeof(company_name)", incorrect_msg="Did you store textual data in `company_name`?")

test_error()
success_msg("Awesome!")

# Other functions to note:
#     - test_or(a,b) -- checks if either test a or test b pass
#     - test_ggplot() -- can check if plots are correct
#     - test_function() -- can also check included parameters
#     - test_loop() -- checking for and while loops
#     - test_library_function('package', not_called_msg='',incorrect_msg='')
#     - test_if_else() -- checking if statements
#     - test_expression_error() -- can check if functions are properly defined
#     - test_operator('operator',), not_called_msg='',incorrect_msg='')
#     - test_function_definition() -- rigorously check defined function
#     - test_data_frame() -- check if dataframe [columns] are equivalent
#     - test_function_result, test_expression_result
```