## Agenda

1. Review of Homeworks
2. A human understanding of regression
3. Dinner break
5. Preprocessing and BoxCox
6. The $K$NN algorithm and the Confusion Matrix

# Homework

## HW1

- We need to work on writing quality.
- We need to work on RMSE intepretation.
- We need to work on using `summary` responsibly.
- We need to work on applying lecture topics to leading questions.
- We would benefit from use of the [`embed-resources`](https://quarto.org/docs/reference/formats/html.html) option in Quarto.

## [HW1 Sols Posted](../hws/wine_of_pnw_sol.html)

<iframe src="../hws/wine_of_pnw_sol.html" height="500" width="800" title="hw1 sols"></iframe> 

## Throwback ThMonday

- I took my old grad ML repo down, but I've restored it.
- [Here](https://github.com/cd-public/cdml15)
- Takeaways:
  - Everything is typeset.
    - Mathematics differentiated from
    - Code block differentiated from
    - Technical writing.
  - No missing assets (e.g. images)
  - Printable.

## HW2

- Think
- Pair
- Share

# Today

## Setup

```{r setup}
library(tidyverse)
library(moderndive)
library(class)
library(caret)
library(fastDummies)
library(thematic)
theme_set(theme_dark())
thematic_rmd(bg = "#111", fg = "#eee", accent = "#eee")
wine <- readRDS(gzcon(url("https://cd-public.github.io/D505/dat/wine.rds"))) %>%
  mutate(lprice = log(price))
```


# Reporting Impact from Regressions

## Correlation

[*Credit:* Modern Dive](https://moderndive.com/5-regression.html)

[http://guessthecorrelation.com/](http://guessthecorrelation.com/) ...my high score is 72 (Jameson 122)

![](./images/correlation.png)

## Calculating correlation

```{r}
wine %>% 
  summarise(cor_p=cor(price,points),cor_lp=cor(lprice,points))
```

## Exercise

1. Calculate the correlation between $\log$(price) and points... 
2. ...by variety...
3. ...for Oregon Chardonnay, Pinot Noir and Pinot Gris...
4. ...in the same tibble!

## Solution

```{r}
wine %>% 
  filter(province=="Oregon") %>% 
  filter(variety %in% c("Chardonnay","Pinot Noir","Pinot Gris")) %>% 
  group_by(variety) %>% 
  summarise(correlation=cor(lprice,points))
```

## Visualizing these different correlations

```{r}
#| output-location: slide
wine %>% 
  filter(province=="Oregon") %>% 
  filter(variety %in% c("Chardonnay","Pinot Noir","Pinot Gris")) %>% 
  ggplot(aes(points,lprice, color=variety)) +
    geom_point(alpha=0.3)+
    facet_wrap(~variety)+
    geom_smooth(method = lm)
    
```


## Graphing residuals (bad)

```{r}
model <- lm(price~points, filter(wine,province=="Oregon"))
get_regression_points(model) %>% 
  ggplot(aes(points, residual))+
    geom_point()
```

## Annotate

![](images/spread.png)

## Graphing residuals (good)

```{r}
model <- lm(lprice~points, filter(wine,province=="Oregon"))
get_regression_points(model) %>% 
  ggplot(aes(points, residual))+
    geom_point()
```


## Interpreting the coefficients

```{r}
model <- lm(lprice~points, filter(wine,province=="Oregon"))
pct = (exp(coef(model)["points"]) - 1) * 100
c(coef(model)["points"],pct) 
```

- We logged the *dependent* variable (price)
  - A 1 point ratings increase = ``r round(pct,2)``\%
  - That is, a percent change in rating to an *absolute* change in the dependent variable. 
- $(e^x - 1)*100$

## Interpreting the coefficients

```{r}
m_yr <- lm(lprice~year, filter(wine,province=="Oregon"))
yr = (exp(coef(m_yr)["year"]) - 1) * 100
c(coef(m_yr)["year"],yr) 
```

- This is a *de facto* measure of inflation.

## Some Examples

::::{.columns}

:::{.column width=50%}

```{r}
for (x in 1:10) {
    print(c(x/100, (exp(x/100)-1) * 100))
}
```

:::


:::{.column width=50%}

```{r}
for (x in 1:10) {
    print(c(x/010, (exp(x/010)-1) * 100))
}
```

:::

::::

## Pretty Print

```{r}
for(v in c("Chardonnay", "Pinot Gris","Pinot Noir")){
  m <- lm(lprice~points, filter(wine,province=="Oregon", variety==v))
  pct <- round((exp(coef(m)["points"]) - 1) * 100,2)
  print(str_c("For ",v,", a 1 point ratings increase leads to a ",pct,"% increase in price."))
}
```

## Summary

- Only the dependent/response variable is log-transformed. 
  - Exponentiate the coefficient.
  - Subtract one from this number
  - Multiply by 100.
- This gives the percent increase (or decrease).

## $\log$ feature

```{r}
model <- lm(price~lpoints, filter(wine,province=="Oregon") %>% mutate(lpoints=log(points)))
model
```

- What does the *sign* (positive or negative) tell us?
- Was $\log$ appropriate here?

## Percentages

```{r}
coef(model)["lpoints"]/100
```

- Since we logged the IV (feature), a 1% ratings increase is a ~3.24 increase in price on average. 
- What are the units on that?

Note: 
$$ 
x/100 
$$

## LogLog (also elasticity)

```{r}
model <- lm(lprice~lpoints, filter(wine,province=="Oregon") %>% mutate(lpoints=log(points)))
model
```

...a 1\% increase in ratings equals a ``r round(coef(model)["lpoints"],2)``\% increase in price on average


## Units

- Change is per one-unit increase in the independent variable. 
  - Here, independent is points.
  - Dependent is price.

## Example

- For every 1% increase in the independent variable...
  - Basically, one point
  - Our dependent variable increases by about 8.3%.
- A $30 bottle of wine scoring 90 would be worth $32.50 as a 91.
```{r}
30 * (1 + 8.3/100)
```

## Graphing points by variety
```{r}
wine %>% 
  filter(province=="Oregon") %>% 
  filter(variety %in% c("Chardonnay","Pinot Noir","Pinot Gris")) %>% 
  ggplot(aes(variety,points))+
    geom_boxplot()
```

## Summary
```{r}
(tmp <- wine %>% 
  filter(province=="Oregon") %>% 
  filter(variety %in% c("Chardonnay","Pinot Noir","Pinot Gris")) %>% 
  group_by(variety) %>% 
  summarise(mean=mean(points)))
```

- What are the percentage differences here?

## Regression
```{r}
model <- lm(points~variety, 
            filter(wine,province=="Oregon",variety %in% c("Chardonnay","Pinot Noir","Pinot Gris")))
get_regression_table(model)
```

- What types of variables are we considering here?


## Assumptions of linear regression

1. **L**inearity of relationship between variables
2. **I**ndependence of the residuals
3. **N**ormality of the residuals
4. **E**quality of variance of the residuals

## Linearity of relationship
[*Credit:* Modern Dive](https://moderndive.com/5-regression.html)

![](./images/non-linear.png)

What would the residuals look like?

## Independence

- Given our original model of 
$$
\log(\text{price})=m*\text{Points}+b
$$ 

- are there any problems with independence?
- How could we check?

## Normality
```{r}
model <- lm(lprice~points, filter(wine,province=="Oregon"))
get_regression_points(model) %>% 
  ggplot(aes(residual))+
    geom_histogram(color="white")
```

## Equality of variance

```{r}
get_regression_points(model) %>% 
  ggplot(aes(points, residual))+
    geom_jitter(alpha=0.2)
```

## No equality in the variance

[*Credit:* Modern Dive](https://moderndive.com/5-regression.html) (click)

![](./images/unequal-variance.png)


# Dinner

![](./images/comic3.png)

# Preprocessing and BoxCox

## Setup

- Pivot to pinot.
  - Say "$\pi^0$"

```{r}
wine <- readRDS(gzcon(url("https://cd-public.github.io/D505/dat/pinot_orig.rds")))
```

## Preprocessing

- Box-Cox transformations use maximum likelihood estimate to estimate value for $\lambda$

$$
y(\lambda) = \frac{x^{\lambda}-1}{\lambda}
$$

- The goal is to make data seem more like a normal distribution.

## in R

- LaTeX

$$
y(\lambda) = \frac{x^{\lambda}-1}{\lambda}
$$

- R
```{r}
y <- function(A) {
  r <- function(x) {
    (x ** A - 1) / A
  }
}
```

## Examples

- When $\lambda=1$, there is no transformation

$$
y(1) = \frac{x^{\lambda}-1}{\lambda} = 
      \frac{x^{1}-1}{1} =
        x-1 \approx x
$$

$$
f = y(-1) \implies f(x) \approx x
$$

## Examples

- When $\lambda=0$, it is log transformed

$$
y(0) = \frac{x^{\lambda}-1}{\lambda} = 
      \frac{x^{0}-1}{0}
$$

$$
f = y(0) \implies f(x) \approx \log(x)
$$


- Zero is a special case, but using a little thing called "calculus" this sorta makes sense.
- Basically, negative infinity at 0, then increase slowly once positive.

## $\lambda = 0$

```{r}
for (x in 1:3) {
  print(x/10)
  for (l in 3:1) {
    f = y(l/100)
    print(c(l/100, x/10, f(x/10), log(x/10), f(x/10) - log(x/10)))
  }
}
```

## Examples

$$
y(.5) = \frac{x^{\lambda}-1}{\lambda} = 
      \frac{x^{.5}-1}{.5} =
      2\times(\sqrt{x}-1) \approx \sqrt{x}
$$

$$
f = y(.5) \implies f(x) \approx \sqrt{x}
$$


## Examples

- When $\lambda=-1$, it is an inverse

$$
y(-1) = \frac{x^{\lambda}-1}{\lambda} = 
      \frac{x^{-1}-1}{-1} =
      \frac{x^{-1}}{-1}+\frac{-1}{-1} =
      -x^{-1}-1 \approx x^{-1}
$$
$$
f = y(-1) \implies f(x) \approx x^{-1}
$$

## A Note

- I am only aware of the following Box-Cox formulation:

$$
y(\lambda) =
\begin{cases}
 \dfrac{y_i^\lambda - 1}{\lambda} & \text{if } \lambda \neq 0, \\
 \ln y_i & \text{if } \lambda = 0,
\end{cases}
$$

- Box-Cox is specified this way on Wikipedia, in Scikit documentation, in Scipy documentation.
- I can only find one reference to $y(1)(x) \equiv x^{-1}$ [here](https://www.statisticshowto.com/probability-and-statistics/normal-distributions/box-cox-transformation/)
  - It's citation are closed source.
- There is no mention in the companion text [ILSR](https://www.statlearning.com/)

## My Theory

::::{.columns}


:::{.column width=50%}

- This chart showed up in literature somewhere
- Either a miscalculation or just some other transform.
- It may be used in practice, including possibly in R?
- It *doesn't matter* they only change wrt scaled values.

:::


:::{.column width=50%}

<table>
<tr>
<td>-3</td>
<td>Y<sup>-3</sup> = 1/Y<sup>3</sup></td>
</tr>
<tr>
<td>-2</td>
<td>Y<sup>-2</sup> = 1/Y<sup>2</sup></td>
</tr>
<tr>
<td>-1</td>
<td>Y<sup>-1</sup> = 1/Y<sup>1</sup></td>
</tr>
<tr>
<td>-0.5</td>
<td>Y<sup>-0.5</sup> = 1/(√(Y))</td>
</tr>
<tr>
<td>0</td>
<td>log(Y)**</td>
</tr>
<tr>
<td>0.5</td>
<td>Y<sup>0.5</sup> = √(Y)</td>
</tr>
<tr>
<td>1</td>
<td>Y<sup>1 </sup>= Y</td>
</tr>
<tr>
<td>2</td>
<td>Y<sup>2</sup></td>
</tr>
<tr>
<td>3</td>
<td>Y<sup>3</sup></td>
</tr></table>

:::

::::

## Just use the function

- R: envstats
```{.r}
library(envstats)

# Example data
data <- c(1, 2, 3, 4, 5)

# Apply Box-Cox transformation
boxcox_result <- boxcoxTransform(data)

# Print the result
print(boxcox_result)
```


## Just use the function

- Py: scipy.stats
```{.py}
from scipy.stats import boxcox
import numpy as np

# Example data
data = np.array([1, 2, 3, 4, 5])

# Apply Box-Cox transformation
boxcox_result, lambda_value = boxcox(data)

# Print the result
print("Transformed data:", boxcox_result)
print("Lambda value:", lambda_value)
```

## On Python

- I like the Python [boxcox documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.boxcox.html):

```{.haskell}
y = (x**lmbda - 1) / lmbda,  for lmbda != 0
    log(x),                  for lmbda = 0
```

- This was how I tracked down what I believe to be the inconsistency with other Box-Cox definitions.

# Onward

## Caret preprocessing is so easy!

```{r}
wine %>% 
  preProcess(method = c("BoxCox","center","scale")) %>% 
  predict(wine) %>% 
  select(-description) %>% 
  head()
```


## Or is it?

But wait... what is wrong here?

```{r}
wino <- wine %>%
  mutate(year_f = as.factor(year))

wino <- wino %>% 
  preProcess(method = c("BoxCox","center","scale")) %>% 
  predict(wino)

head(wino %>% select(starts_with("year")))
```

- Are years normally distributed?


# The $K$NN Algorithm

## Algorithm {.smaller}

1. Load the data
2. Initialize $K$ to your chosen number of neighbors
3. For each example in the data
  - Calculate the distance between the query example and the current example from the data.
  - Add the distance and the index of the example to an ordered collection
4. Sort the ordered collection of distances and indices from smallest to largest (in ascending order) by the distances
5. Pick the first $K$ entries from the sorted collection
6. Get the labels of the selected $K$ entries
7. If regression, return the mean of the $K$ labels
8. If classification, return the mode of the $K$ labels

## Basis

- We assume:
  - Existing datapoints in something we think of as a space
    - That is, probably two numerical value per point in a coordinate plane
    - Categorical is fine - think a Punnett square
  - Existing datapoints are labelled
    - Numerical or categorical still fine!
- To visualize, we will have a 2d space with color labels.

## Let's draw it

<a  title="Antti Ajanki AnAj, CC BY-SA 3.0 &lt;http://creativecommons.org/licenses/by-sa/3.0/&gt;, via Wikimedia Commons" href="https://commons.wikimedia.org/wiki/File:KnnClassification.svg"><img style="background-color:white;" width="50%" alt="KnnClassification" src="https://upload.wikimedia.org/wikipedia/commons/thumb/e/e7/KnnClassification.svg/256px-KnnClassification.svg.png?20101222011435"></a>

## Let's draw it 2


<a title="Paolo Bonfini, 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:KNN_decision_surface_animation.gif"><img width="100%" alt="KNN decision surface animation" src="https://upload.wikimedia.org/wikipedia/commons/thumb/7/78/KNN_decision_surface_animation.gif/512px-KNN_decision_surface_animation.gif?20240718130312"></a>

## Let's draw it 3


![](https://media.geeksforgeeks.org/wp-content/uploads/20231207103856/KNN-Algorithm-(1).png)

[G4G](https://www.geeksforgeeks.org/k-nearest-neighbours/)

## Engineering some features

- Create an "other" for most tasters.

```{r}
wino <- wino %>% 
  mutate(taster_name = fct_lump(taster_name,5))
```

## Engineering some features

- Create dummys for years, tasters
```{r}
wino <- wino %>% dummy_cols(
    select_columns = c("year_f","taster_name"),
    remove_most_frequent_dummy = T, 
    remove_selected_columns = T)
```

## Engineering some features

- Convert everything to snake case.

```{r}
wino <- wino %>% 
  rename_all(funs(tolower(.))) %>% 
  rename_all(funs(str_replace_all(., "-", "_"))) %>% 
  rename_all(funs(str_replace_all(., " ", "_")))
```

## Engineering some features

- Add indicators for 3 tasting notes.

```{r}
wino <- wino %>% 
  mutate(note_cherry = str_detect(description,"cherry")) %>% 
  mutate(note_chocolate = str_detect(description,"chocolate")) %>%
  mutate(note_earth = str_detect(description,"earth")) %>%
  select(-description)
```

## Engineering some features

- Let's see it

```{r}
head(wino) %>% 
  select(1:6)
```

## Split

```{r}
set.seed(505)
wine_index <- createDataPartition(wino$province, p = 0.8, list = FALSE)
train <- wino[ wine_index, ]
test <- wino[-wine_index, ]
```

## Simple model

- Specify a $K$NN model.

```{r}
fit <- knn(
  train = select(train,-province), 
  test = select(test,-province), 
  k=5, 
  cl = train$province)
```

## Confusion matrix

- Let's look at Kappa.

```{r}
confusionMatrix(fit,factor(test$province))$overall
```

## [Kappa $\kappa$ statistic](https://pmc.ncbi.nlm.nih.gov/articles/PMC4916348/)

> Kappa statistic is a measurement of the agreement for categorical items Kappa can be used to assess the performance of kNN algorithm.

$$
\kappa = \dfrac{P(A)-P(E)}{1 - P(E)}
$$

> where $P(A)$ is the relative observed agreement among raters, and $P(E)$ is the proportion of agreement expected between the classifier and the ground truth by chance.

## Kappa $\kappa$ statistic

Rule of thumb.

- \< 0.2 (not so good)
- 0.21 - 0.4 (ok)
- 0.41 - 0.6 (pretty good)
- 0.6 - 0.8 (great)
- \> 0.8 (almost perfect)

We had ~.9...


# Overfitting... or a leak?

## Review the dataframe

```{r}
summary(wino)
```

## Determine what dominates

```{r}
omit <- function(prefix) {
  train <- select(train, -starts_with(prefix))
  test <- select(test, -starts_with(prefix))
  fit <- knn(
    train = select(train,-province), 
    test = select(test,-province), 
    k=5, 
    cl = train$province)
  confusionMatrix(fit,factor(test$province))$overall
}
```

## Test

```{r}
c(omit('year')["Kappa"], omit('note')["Kappa"], omit('taster')["Kappa"])
```

## Fixing the leak

- Dastardly humans, always existing in a physical location.

```{r}
train <- select(train, -starts_with("taster"))
test <- select(test, -starts_with("taster"))
```

## Rerun

```{r}
fit <- knn(
  train = select(train,-province), 
  test = select(test,-province), 
  k=5, 
  cl = train$province)
```

- We should probably have written function here!
- That is a lot of lines to copy+paste...

## Confusion matrix

```{r}
confusionMatrix(fit,factor(test$province))
```

## With parameter tuning over $K$

```{r}
fit <- train(province ~ .,
             data = train, 
             method = "knn",
             tuneLength = 15,
             trControl = trainControl(number = 1)) # default bootstrap
fit
```


## Confusion Matrix
```{r}
confusionMatrix(predict(fit, test),factor(test$province))
```


## Tuning and subsampling

```{r}
fit <- train(province ~ .,
             data = train, 
             method = "knn",
             tuneLength = 15,
             metric = "Kappa", # this is new
             trControl = trainControl(number = 1))
fit
```

## Tuning plot

```{r}
ggplot(fit, metric="Kappa")
```


## Group modeling problem I

* Practice running different versions of the model
* Create some new features and...
* See if you can achieve a Kappa >= 0.5!

$$
\kappa \geq 0.5
$$


## Bonus: KNN for regression

```{r}
fit <- train(price ~ .,
             data = train, 
             method = "knn",
             tuneLength = 15,
             trControl = trainControl(number = 1))
fit

```

