## Agenda

## Agenda

0. Course Announcements
1. Mathematics of logistic regression
2. Implementation with Caret
4. ROC Curves
5. Group work


# Midterm 3/17

## Brief recap of course so far {.smaller}

- Linear regression (e.g. `lprice ~ .`), assumptions of model, interpretation.
- $K$-NN (e.g., `province ~ .`), multi-class supervised classification. Hyperparameter $k$.
- Naive Bayes (e.g., `province ~ .`), multi-class supervised classification.
- Logistic regression (e.g., `province=="Oregon" ~ .`), binary supervised classification. Elastic net. 
- Feature engineering (logarithms, center/scaling, Box Cox, tidytext, etc.).
- Feature selection (correlation, linear / logistic coefficients, frequent words, frequent words by class, etc.).

## Practice

- I am working to develop a practice midterm.
- I will circulate it by 3 Mar.
- You will go over it on 10 Mar.
- It is based on the 5 homeworks.
- It is based on the prior slide.
  - Little to no computatational linguistics
  - I'm regarding `tidytext` as extension, not core, content.

## Modality Discussion

- I would release an assignment electronically at 6 PM
- We can do in person or otherwise.
- It will be "cheat proof"
  - I will ask you nothing for which it matters how you determine the answer.
  - If e.g. ChatGPT can be mind controlled into doing high quality feature engineering, you get points for mind controlling ChatGPT.

# First Model Due 3/10

##  Publish

Each group should create:
1. An annotated `.*md` file, and 
2. The .rds/.pickle/.parquet file that it generates, that
3. Contains *only* the features you want in the model.

Under version control, on GitHub.

##  Constraints

I will run:

1. The specified $K$NN or Naive Bayes model,
2. With: `province ~ .` (or the whole data frame in `scikit`)
3. With repeated 5-fold cross validation
4. With the same index for partitioning training and test sets for every group.
5. On whatever is turned in before class.
6. Bragging rights for highest Kappa

## Context

- The "final exam" is that during the last class you will present your model results as though you are speaking to the
managers of a large winery.
- It should be presented from a Quarto presentation on GitHub or perhaps e.g. RPubs.
- You must present via the in-room "teaching machine" computer, not your own physical device, to ensure that you are comfortable distributing your findings.

## Group Meetings

- You should have a group assignment
- Meet in your groups!
- Talk about your homework *with* your group.

# Logistic Regression

## Algorithm

- Assume a linear relationship between the log odds and a set of predictor variables.

$$
log(\frac{p}{1-p})=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}
$$

With a bit of algebra you can get the probabilities as...

$$
p=\frac{1}{1+e^{-(\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2})}}
$$

## 
$$
log(\frac{p}{1-p})=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}
$$

'p' represents the probability of the event occurring.

- Represented as natural log `ln`...
- of the odds ratio ($\frac{p}{1-p}$)...
- as a linear combination of 
  - predictor variables ($x_1$, $x_2$) and
  - their corresponding coefficients ($\beta_1$, $\beta_2$) 
  - plus an intercept ($\beta_0$). 

## Understanding check?

- Why do we call this regression instead of classification?
- Quick Algebra review
  - Look, it's better to be comfortable with algebra than not.

##

<a title="Canley, CC BY-SA 4.0 &lt;https://creativecommons.org/licenses/by-sa/4.0&gt;, via Wikimedia Commons" href="https://commons.wikimedia.org/wiki/File:Exam_pass_logistic_curve.svg"><img width="100%" alt="Logistic regression curve showing probability of passing an exam versus hours studying" src="https://upload.wikimedia.org/wikipedia/commons/thumb/c/cb/Exam_pass_logistic_curve.svg/512px-Exam_pass_logistic_curve.svg.png?20220327133556"></a>


## Log Odds

- **1. Start with the log-odds equation:**

$$
log(\frac{p}{1-p}) = \beta_0 + \beta_1x_1 + \beta_2x_2
$$



## Exponentiate

- **2. Exponentiate both sides:**

- To remove the logarithm, we take the exponential (e) of both sides of the equation:

$$
e^{log(\frac{p}{1-p})} = e^{\beta_0 + \beta_1x_1 + \beta_2x_2}
$$

## Simplify

- **3. Simplify the left side:**

- The exponential function and the natural logarithm are inverse functions of each other.  
- Therefore, $e^{log(x)} = x$.  This simplifies the left side:

$$
\frac{p}{1-p} = e^{\beta_0 + \beta_1x_1 + \beta_2x_2}
$$

## Rewrite

- **4. Rewrite the right side using exponent rules:**

- We can rewrite the right side using the rule $e^{a+b} = e^a * e^b$:

$$
\frac{p}{1-p} = e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2}
$$

## Find $p$ 1

* 5. Find $p$
  * **Multiply both sides by (1-p):**

$$
p = (1-p) * e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2}
$$

## Find $p$ 2

* 5. Find $p$
  * **Distribute the right-hand side (RHS):**

$$
p = e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2} - p * e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2}
$$

## Find $p$ 2

* 5. Find $p$
  * **Move the $p$ term to the left side:**

$$
p + p * e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2} = e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2}
$$

## Find $p$ 3

* 5. Find $p$
  * **Factor out $p$ on the left side:**

$$
p * (1 + e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2}) = e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2}
$$

## Find $p$ 4

* 5. Find $p$
  * **Divide both sides by the term in parentheses:**
$$
p = \frac{e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2}}{1 + e^{\beta_0} * e^{\beta_1x_1} * e^{\beta_2x_2}}
$$

## Simplify

- **6. Simplify the expression:**

- We can further simplify by multiplying the numerator and denominator by $e^{-(\beta_0 + \beta_1x_1 + \beta_2x_2)}$:

$$
p = \frac{1}{e^{-(\beta_0 + \beta_1x_1 + \beta_2x_2)} + 1}
$$

## Rewrite

- Which is commonly written as:

$$
p = \frac{1}{1 + e^{-(\beta_0 + \beta_1x_1 + \beta_2x_2)}}
$$


# Implementation with Caret

## Libraries Setup
```{r libs}
sh <- suppressPackageStartupMessages
sh(library(tidyverse))
sh(library(caret))
sh(library(tidytext))
sh(library(SnowballC)) # new?
sh(library(pROC))      # new?
sh(library(glmnet))
data(stop_words)
sh(library(thematic))
theme_set(theme_dark())
thematic_rmd(bg = "#111", fg = "#eee", accent = "#eee")
```

## Dataframe

```{r fg}
wine <- readRDS(gzcon(url("https://cd-public.github.io/D505/dat/pinot.rds")))
names(wine)[names(wine) == 'id'] = 'id'
```


## Extract certain words

- We'll proceed in a few steps.
  - First, a function to extract all words.
    - We'll omit stopwords
      - Remember those? Like "of" and "a"
    - We'll omit a few words-of-choice
      - Like "pinot"

## Extract all words

```{r desc_to_words}
desc_to_words <- function(df, omits) { 
  df %>%
    unnest_tokens(word, description) %>%
    anti_join(stop_words) %>% # get rid of stop words
    filter(!(word %in% omits))
}
```

## Extract all  words

```{r run desc_to_words}
words <- desc_to_words(wine, c("wine","pinot","vineyard"))
# The second argument is our custom stopwards, as a vector
head(words)
```

## Extract certain words

- We'll proceed in a few steps.
  - Optionally, look at stems
    - This is really cool.
    - Do I care about the difference between "acidic" and "acidity"?
      - The linguistics consensus leans toward no.

## STEM

- Short for data *S*cience, da*T*a science, data scienc*E*, and *M*achine learning.
```{r words_to_stems}
words_to_stems <- function(df) { 
  df %>%
    mutate(word = wordStem(word))
}
```
- We want to use a `tidytext` built-in here, of course.
  - We aren't domain experts in linguistics...
  - But we are domain experts in using R libraries.

## Example

> [Stem may either consist of a root (e.g. run) alone or a compound word, such as meatball and bottleneck (examples of compound nouns) or blacken and standardize (examples of compound verbs). The stem of the verb to wait is wait: it is the part that is common to all its inflected variants.](https://en.wikipedia.org/wiki/Word_stem#Root_vs_stem)

## Wait

- 'wait' (infinitive, imperative, present subjunctive, and present indicative except in the 3rd-person singular)
- 'waits' (3rd person singular simple present indicative)
- 'waited' (simple past)
- 'waited' (past participle)
- 'waiting' (present participle)

## STEM

```{r run words_to_stems}
stems <- words_to_stems(words)
head(stems)
```

## Bottl

- Another example.

```{r g6}
g6_df <- tibble(description = "Popping bottles in the ice; like a blizzard.")
words_to_stems(desc_to_words(g6_df, c("ice")))
```

## Like a G6

<iframe width="100%" height="100%" src="https://www.youtube.com/embed/w4s6H4ku6ZY?si=DWmx3skbebkcm1Dk" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture; web-share" referrerpolicy="strict-origin-when-cross-origin" allowfullscreen></iframe>


## Word Count

- With either or words or our stems, we can see how many of each word we have easily enough.
  - And eliminate words will less than a certain count!
```{r filter_by_count}
filter_by_count <- function(df, j) { 
  df %>%
    count(id, word) %>% 
    group_by(id) %>% mutate(exists = (n>0)) %>% ungroup %>% 
    group_by(word) %>% 
    mutate(total = sum(n)) %>% 
    filter(total > j)
}
```

## Check it 

:::: {.columns}

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

```{r check words}
fwords = filter_by_count(words, 1000)
head(fwords)
```

:::

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

```{r check stems}
fstems = filter_by_count(stems, 1000)
head(fstems)
```

:::

::::

## Back to Wine

- We currently have multiple entries per ID in a tidy dataframe
- We would like to get back to have features for each ID, as each ID is some wine.
- One way to do so, is to make columns from the word data.
- We use `pivot_wider`

## Pivot

```{r pivoter}
pivoter <- function(words, df) {
  words %>%
    pivot_wider(id_cols = id, names_from = word, values_from = exists, values_fill = list(exists=0)) %>% 
    right_join(select(df,id,province)) %>% 
    drop_na() %>% 
    select(-id)
}
```

## Check it 

:::: {.columns}

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

```{r pivot words}
head(pivoter(fwords, wine))
```

:::

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

```{r pivot stems}
head(pivoter(fstems, wine))
```

:::

::::

## Full Function

- Create a function to extract words with totals > j

```{r wordfunc}
wine_words <- function(df, j, stem) { 

  words <- desc_to_words(df, c("wine","pinot","vineyard"))
  
  if (stem) {
    words <- words_to_stems(words)
  }
  
  words <- filter_by_count(words, j)

  pivoter(words, df)
}
```




## Look at the data

```{r w100}
wino <- wine_words(wine, 1000, F)

wino %>% 
  head(10) %>% 
  select(1:5, province)
```

## Logistic Regression

- We check a true/false outcome.

```{r w101}
wino <- wino %>% 
  mutate(oregon = factor(province=="Oregon")) %>%
  select(-province)

wino %>% 
  head(10) %>% 
  select(1:5, oregon)
```

## Split the data 

```{r tt100}
wine_index <- createDataPartition(wino$oregon, p = 0.80, list = FALSE)
train <- wino[ wine_index, ]
test <- wino[-wine_index, ]
table(train$oregon)
```

## A basic model

```{r getfit}
control = trainControl(method = "cv", number = 5)
get_fit <- function(df) {
  train(oregon ~ .,
        data = df, 
        trControl = control,
        method = "glm",
        family = "binomial",
        maxit = 5) # speed it up - default 100
}
fit <- get_fit(train)
```

## Check Kappa

```{r}
fit
```

## Probability

- See top coefficients

```{r getodds}
get_odds <- function(fit) {
  as.data.frame(t(exp(coef(fit$finalModel))))   %>%
  rownames_to_column(var = "name") %>%
  pivot_longer(-name, names_to = "class", values_to = "odds") %>%
  arrange(desc(odds)) %>%
  head()
}
get_odds(fit)
```

## Confusion Matrix

```{r getmatrix}
get_matrix <- function(fit, df) {
  pred <- factor(predict(fit, newdata = df))
  confusionMatrix(pred,factor(df$oregon))
}
get_matrix(fit,test)
```

Not bad. But what if we decrease threshold for a word to be included?

## Using more words

```{r w500}
wino <- wine_words(wine, 500, F) %>% 
  mutate(oregon = factor(province=="Oregon")) %>%
  select(-province)
wine_index <- createDataPartition(wino$oregon, p = 0.80, list = FALSE)
train <- wino[ wine_index, ]
test <- wino[-wine_index, ]

fit <- get_fit(train)
```

## Check Kappa

```{r}
fit
```

## Check odds

```{r}
get_odds(fit)
```

## Confusion Matrix

```{r 500matrix}
get_matrix(fit,test)
```

## Using stems

```{r stem}
wino <- wine_words(wine, 1000, T) %>% 
  mutate(oregon = factor(province=="Oregon")) %>%
  select(-province)
wine_index <- createDataPartition(wino$oregon, p = 0.80, list = FALSE)
train <- wino[ wine_index, ]
test <- wino[-wine_index, ]

fit <- get_fit(train)
```

## Check Kappa

```{r}
fit
```

## Check odds

```{r}
get_odds(fit)
```

## Confusion Matrix

```{r stematrix}
get_matrix(fit,test)
```

Even better?

# Receiver Operating Characteristic (ROC) Curve

## ROC Curve

> [A receiver operating characteristic curve, or ROC curve, is a graphical plot that illustrates the performance of a binary classifier model (can be used for multi class classification as well) at varying threshold values.](https://en.wikipedia.org/wiki/Receiver_operating_characteristic)

## Image

![](https://upload.wikimedia.org/wikipedia/commons/6/6b/Roccurves.png)

## Sensitivity

- **Sensitivity**: True Positive Rate 
- Measures how well the model detects true positives.
  - TP: True Positive
  - FN: False Negative
  
$$ 
\text{Sensitivity} = \frac{\text{TP}}{\text{TP} + \text{FN}}
$$ 

## Specificifity

- **Specificity**: True Negative Rate
- Measures how well the model avoids false positives.
  - FP: False Positive
  - TN: True Negative

$$  
\text{Specificity} = \frac{\text{TN}}{\text{TN} + \text{FP}}
$$  

## Storm Prediction

- **True Positive (TP):** The model predicts a storm, and a storm actually occurs.
- **False Positive (FP):** The model predicts a storm, but no storm occurs (false alarm).
- **True Negative (TN):** The model predicts no storm, and no storm occurs.
- **False Negative (FN):** The model predicts no storm, but a storm actually occurs (missed event).

## Goals

- A perfect model would have both 
  - high sensitivity 
    - (detecting all real storms) and 
  - high specificity 
    - (avoiding false alarms). 
- The ROC curve helps analyze this trade-off.

## Interpreting the ROC Curve

- A point near (0,1) on the ROC curve represents high sensitivity and high specificity (ideal performance).

![](https://upload.wikimedia.org/wikipedia/commons/3/36/ROC_space-2.png)

## Interpreting the ROC Curve

- A curve closer to the **diagonal line** (random guessing) indicates poor predictive ability.

![](https://upload.wikimedia.org/wikipedia/commons/3/36/ROC_space-2.png)

## Interpreting the ROC Curve

- The **Area Under the Curve (AUC)** quantifies the overall model performance, with **AUC = 1** being perfect and **AUC = 0.5** being no better than chance.

![](https://upload.wikimedia.org/wikipedia/commons/1/13/Roc_curve.svg)

## Model Tuning

- In e.g. a storm model....
- adjusting the model's decision threshold 
  -  (e.g., how strong a weather signal needs to be before predicting a storm), 
  - we can move along the ROC curve to balance 
    - sensitivity and 
    - specificity.

##  ROC Curve evaluation

- Let's look at our most recent fit:

```{r numeric}
fit
```

## See the curve

```{r}
prob <- predict(fit, newdata = test, type = "prob")[,2]
myRoc <- roc(test$oregon, prob)
plot(myRoc)
auc(myRoc)
```

## Exercise

1. Gather into your prediction teams.
2. Choose a Pinot province other than Oregon or California
3. Use logistic regression to find the words/terms that increase the odds of choosing that province the most

# Vocabulary

- Logistic Regression
- Odds
- Sensitivity
- Specificity
- ROC
- AUC