## Agenda

1.  Review of Homework 1
2.  Feature Engineering I
3.  Dinner Break
4.  The Caret framework
5.  Vocabulary

## Packages

-   Today I use the following libraries:

``` r
local({r <- getOption("repos")
       r["CRAN"] <- "https://cran.r-project.org" 
       options(repos=r)
})
# Old
# install.packages("tidyverse")
# install.packages("caret")
# New?
install.packages("fastDummies")
# Just for the slides
# install.packages("thematic")
```

-   You will have some but perhaps not others.

## Libraries

-   I'll just include them upfront.

```{r Libraries}
library(tidyverse)
library(caret)
library(fastDummies)
# Just for the slides
library(thematic)
theme_set(theme_dark())
thematic_rmd(bg = "#111", fg = "#eee", accent = "#eee")
```

## Setup

-   We will work with a `wine` dataset that is enormous.
    -   Just to render a bit quickly, take a sample.
    -   You are welcome to work with the [full dataset](https://cd-public.github.io/courses/rmls25/dat/wine.rds)!

```{r Setup}
wine <- readRDS(gzcon(url("https://cd-public.github.io/D505/dat/wine.rds")))
# performance concession
# wall <- wine
# wine = wall[sample(nrow(wall), 100), ]
summary(wine)
```

# Feature Engineering

## Group Exercise (30m)

1.  Identify 3 "interesting" features of the wine dataset
2.  **Bonus** 
    -   &exist; Critic "Roger Voss"
    -   &exist; a wine varietal(s) "Voss" seems to dislike
    -   Find said varietal.

## Categorical vs. Continuous

-   What is a categorical variable?
-   What is a continuous variable?
-   Why visualize at the data before modeling it?

## Categorical Example 1

```{r Cat1}
wine %>%
  mutate(roger = taster_name == "Roger Voss") %>%
  mutate(pinot_gris = variety == "Pinot Gris") %>%
  drop_na(roger) %>%
  group_by(roger, pinot_gris) %>%
  summarize(points = mean(points)) %>%
  ggplot() +
  aes(pinot_gris, points, color = roger) +
  geom_line(aes(group = roger), size = 2)
```

## Categorical Example 2

```{r Cat2}
wine %>%
  filter(province == "Oregon") %>%
  group_by(year) %>%
  summarize(price = mean(price)) %>%
  ggplot(aes(year, price)) +
  geom_smooth() +
  labs(title = "Oregon wine over the years")
```

## Exercise (15 min)

1.  Group by winery and year, Find:
    -   The average score, and 
    -   Number of reviews.
2.  Find year-on-year change in score by winery.
3.  How might you use this in prediction? 
    -   What kind of problem might it help with?

## Year-on-Year Change Example

```{r YoY}
wine %>%
  group_by(winery, year) %>%
  summarize(avg_score = mean(points), num_reviews = n_distinct(id)) %>%
  select(year, winery, num_reviews, avg_score) %>%
  arrange(winery, year) %>%
  mutate(score_change = avg_score - lag(avg_score)) %>%
  drop_na(score_change) %>%
  summarize(mean(score_change))
```

## Dummy Variables

-   **What are Dummy Variables?**: 
    -   Represent categories as 0s and 1s for models.\
-   **Why Use Dummy Variables?**: 
    -   Handle categorical data in numerical algorithms.\
-   **Avoid Dummy Trap**: 
    -   Drop one column to prevent multicollinearity.

## Many vs Few Dummies

-   **Few Dummies**: 
    -   Simplifies models, risks losing fine-grained patterns.
-   **Many Dummies**: 
    -   Captures detailed trends, increases model complexity.
-   **Key Decision**: 
    -   Balance interpretability and predictive power.

## "fastDummies" Package

-   **Purpose**: 
    -   Quickly create dummy variables in R datasets.
-   **Key Functions**: 
    -   `dummy_cols()` adds dummy columns efficiently.
-   **Features**: 
    -   Handles multiple columns and missing data flexibly.

## Few Dummies

```{r Few Dummies}
wine %>%
  select(taster_name) %>%
  dummy_cols() %>% # library(fastDummies)
  select(1:4) %>%
  head()
```

## Many Dummies

```{r Many Dummies}
wine %>%
  select(variety) %>%
  mutate(variety = fct_lump(variety, 4)) %>%
  dummy_cols() %>%
  head()
```

## Other types of engineered categorical features...

-   Words or phrases in text
-   A given time period
-   An arbitrary numerical cut-off
-   Demographic variables

## What about numerical features?

```{r Price gg}
wine %>%
  ggplot(aes(price)) +
  geom_histogram()
```

## Take the natural log

```{r Log Price gg}
wine %>%
  ggplot(aes(log(price))) +
  geom_histogram()
```

## Standardizing

-   Create a common scale across variables.
    -   Mean-centering $x-\bar{x}$
    -   Scaling: $x/std(x)$

-   Helps reduce bias when interactions are included. 
    -   (i.e. eliminates variance inflation).

## Other transformations.

-   I use logs \> 95% of the time, standarizing \~40%.
-   There are [many other transformations](http://www.feat.engineering/numeric-one-to-many.html):
    -   YoY, QoQ, etc. (absolute and percent)
    -   log
    -   polynomial transforms
    -   lags!

## Standardize

-   `list(normalized = ~(scale(.) %>% as.vector))`: :
    -   `scale(.)`: Standardizes the "points" column.
    -   `%>% as.vector`: Converts back to a vector.

```{r Normalize}
wine %>% mutate_at("points", list(standardized = ~ (scale(.) %>% as.vector())))
```

## Interaction effects

[This chapter](http://www.feat.engineering/detecting-interaction-effects.html) has a good overview of interactions.

-   Start with domain knowledge.
-   Use visualizations.
-   3-way interactions exist, but are rare.
    -   If you suspect a 3-way, also suspect your suspicions.
    -   Complexity increases exponentially in "ways".
    -   These are notoriously hard to explain.

## Dinner (and virtual high fives)

![](images/comic2.gif)

# Le Bibliothèques _caret_

## Philosophy

```{dot Philosophy}
//| echo: false
digraph G {
    
    bgcolor="#101010";

    node [
        fontcolor = "#e0e0e0",
        color = "#e0e0e0",
    ]

    edge [
        color = "#e0e0e0",
        fontcolor = "#e0e0e0"
    ]
    node [shape=circle];
    A [label="All Data"];

    node [shape=pentagon];
    B [label="Training"];
    C [label="Testing"];

    node [shape=rectangle];
    D [label="Resample 1"];
    E [label="Resample 2"];
    F [label="Resample B"];

    node [shape=ellipse];
    G [label="Analysis"];
    H [label="Assessment"];
    I [label="Analysis"];
    J [label="Assessment"];
    K [label="Analysis"];
    L [label="Assessment"];

    A -> B;
    A -> C;
    B -> D;
    B -> E;
    B -> F;
    D -> G;
    D -> H;
    E -> I;
    E -> J;
    F -> K;
    F -> L;
}
```

## Types of resampling

-   [V-fold Cross-Validation](http://www.feat.engineering/resampling.html#cv)
    -   Divides data into $k$ folds, trains on $k−1$ folds, validates on the remaining fold, for all folds.
-   [Monte Carlo Cross-Validation](http://www.feat.engineering/resampling.html#monte-carlo-cross-validation)
    -   Randomly splits data into training and validation sets multiple times, averaging results for evaluation.
-   [The Bootstrap](http://www.feat.engineering/resampling.html#the-bootstrap)
    -   Uses resampling with replacement to estimate model accuracy and variability.

## Setup the Dataframe

-   Follow [this link](https://topepo.github.io/caret) for the full documentation on caret.

```{r Engineering}
wino <- wine %>% # 3 engineered features
  mutate(fr = (country == "France")) %>%
  mutate(cab = str_detect(variety, "Cabernet")) %>%
  mutate(lprice = log(price)) %>%
  drop_na(fr, cab) %>%
  select(lprice, points, fr, cab)
```

-   Off hand, I would've standarized points as well, but
-   We're following Jameson's code...
    -   ...who *understands the data better*.

## Split Samples

```{r Split}
wine_index <- createDataPartition(wino$lprice, p = 0.8, list = FALSE)
wino_tr <- wino[wine_index, ]
wino_te <- wino[-wine_index, ]
summary(wino_tr)
```

## Train the model

-   Configure `train` to cross validate

```{r Train Model}
m1 <- train(lprice ~ .,
  data = wino_tr,
  method = "lm",
  trControl = trainControl(method = "repeatedcv", number = 5, repeats = 3)
)
m1
```

## RMSE outputs

```{r RMSE}
print(m1$resample)
```

## Train vs. test

::::: columns

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

```{r Train}
m1
```

:::

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

```{r Test}
postResample(pred = predict(m1, wino_te), obs = wino_te$lprice)
```

:::

:::::

## Group Exercise (30+ minutes)

1.  Create 5-10 new features (in addition to points)
2.  Create training and test data
3.  For each, train a linear model for log(price)
4.  Report RMSE on test set and cross-validated score.
5.  (Re-)Engineer new(ish) features to lower the RMSE.

# Feature selection

## Stepwise selection

-   **What is Stepwise Selection?**: Iterative method to add or remove predictors in a model based on statistical criteria.\
-   **Types**: Forward selection starts with no predictors; backward elimination starts with all predictors; stepwise combines both.\
-   **Goal**: Identify a model with strong predictive power and minimal overfitting.

## Stepwise selection is bad

Harrell (2015) provides a comprehensive indictment of the method that can be encapsulated by the statement:

> **“… if this procedure had just been proposed as a statistical method, it would most likely be rejected because it violates every principle of statistical estimation and hypothesis testing.”**

  Reference: [Harrell, F. 2015. Regression Modeling Strategies. Springer.](https://link.springer.com/book/10.1007/978-3-319-19425-7https://link.springer.com/book/10.1007/978-3-319-19425-7)

## Engineer 9 features

```{r Mutates}
wino <- wine %>%
  mutate(country = fct_lump(country, 4)) %>%    # 1:4,
  mutate(variety = fct_lump(variety, 4)) %>%    # 5:8,
  mutate(lprice = log(price)) %>%               #   9
  select(lprice, points, country, variety) %>%
  drop_na(.)
head(wino)
```

## Add Dummy Columns
- Careful - a destructive update to `wino`!
```{r Format}
renamer <- function(s) {
  s %>% tolower() %>% str_replace("-| ", "_")
}

wino <- wino %>%
  dummy_cols(remove_selected_columns = TRUE) %>%
  rename_with(.fn = renamer) %>%
  select(-ends_with("other"))
head(wino)
```

## Basic Model
-   Partition
```{r Basic}
wine_index <- createDataPartition(wino$lprice, p = 0.8, list = FALSE)
wino_tr <- wino[wine_index, ]
wino_te <- wino[-wine_index, ]
```
- We would model the same way, so let's take aside.

## Aside: Factoring
- Same modelling command
```
mx <- train(lprice ~ .,
  data = wino_tr,
  method = "lm",
  trControl = trainControl(method = "repeatedcv", number = 5, repeats = 3)
)
```
-   I should factor this into a function.
```{r Factor}
do_training <- function(df, formula) {
  train(formula,
    data = df,
    method = "lm",
    trControl = trainControl(method = "repeatedcv", number = 5, repeats = 3)
  )
}
```

## Train vs. test

::::: columns

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

```{r Train 2}
m2 <- do_training(
  wino_tr, lprice ~ .
)
m2
```

:::

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

```{r Test 2}
postResample(
  pred = predict(m2, wino_te),
  obs = wino_te$lprice
)
```
:::

:::::


## Variable Importance
-   Importance depends on model used...

```{r Importance}
plot(varImp(m2, scale = TRUE))
```

## Variable Importance

-   Each (linear model) coefficient has a standard error, 
    -   Measures certainty of coefficient given data.
-   For the t-statistic, 
    -   Confidence that the coefficient is different from 0
    -   We divide the coefficient by the standard error.
-   If "small" error relative to coefficient
    -   Then "big" t-statistic & high feature importance!
-   What about coefficient as variable importance?

## [Recursive Feature Elimination](https://topepo.github.io/caret/recursive-feature-elimination.html) {.smaller}

1. Tune/train the model on the training set using all predictors.
2. Calculate model performance.
3. Calculate variable importance or rankings.
4. **for** each subset size $S_i$, i = 1...S **do**
    1. Keep the $S_i$ most important variables.
    2. [Optional] Pre-process the data.
    3. Tune/train the model on the training set using $S_i$ predictors.
    4. Calculate model performance.
    5. [Optional] Recalculate the rankings for each predictor.
5. **end**
6. Calculate the performance profile over the $S_i$.
7. Determine the appropriate number of predictors.
8. Use the model corresponding to the optimal $S_i$.

## Size Drop
-   It did not seem like 2024 `r` could handle 90k wine samples.
-   We drop to the 1k data set for this demonstration.
-   There are ways to address performance.
    -   I say: Out-of-scope.
```{r 1k}
wino <- wino[sample(nrow(wino), 1000), ]
```

## Partition Again
-   Partition
```{r Repartition}
wine_index <- createDataPartition(wino$lprice, p = 0.8, list = FALSE)
wino_tr <- wino[wine_index, ]
wino_te <- wino[-wine_index, ]
```

## Caret RFE

```{r RFE}
control <- rfeControl(functions = rfFuncs, method = "cv", number = 2)
# run the RFE algorithm
results <- rfe(select(wino_tr, -lprice), wino_tr$lprice, sizes = c(1:3), rfeControl = control)
# summarize the results
print(results)
```

## See Results
```{r See Results}
predictors(results)
```

## Plot Results

```{r Plot Results}
plot(results)
```

## Practical Workflow

```{dot Practice}
//| echo: false
digraph feature_engineering_pipeline {
    
    bgcolor="#101010";

    node [
        fontcolor = "#e0e0e0",
        color = "#e0e0e0",
    ]

    edge [
        color = "#e0e0e0",
        fontcolor = "#e0e0e0"
    ]
    node [shape=box];
    "Raw Data" -> "Lots of Features" [label="Feature Engineering"];
    "Lots of Features" -> "Candidate Features" [label="Feature Selection"];
    "Candidate Features" -> "Shortlist Features" [label="Expert Input"];
    "Shortlist Features" -> "Finalist Models" [label="DS Judgement"];
    "Finalist Models" -> "Production" [label="Business Unit"];
}
```

## Key Terms

::::: columns

::: {.column width="50%"}
-   Feature Engineering
-   Categorical Feature
-   Continuous Feature
-   Dummy
-   Interaction
:::

::: {.column width="50%"}
-   Caret
-   Model
-   Resampling
-   Train vs. Test Data
-   Variable Importance
:::

:::::


# Bonus Slides:<br> Linear Regression

## 5 Assumptions of Linear Regression

-   Linear regressions have a well-developed statistical theory.

-   This brings perks like confidence intervals on predictions.

-   It also has "costs" in that assumptions need to be satisfied.

## The Five

1. Linearity
2. Constant variance
3. Normality
4. Imperfect multicollinearity
5. Exogeneity

## 1.  Linearity 

-   **The dependent variable is a linear combination of the features.**

-   This is less of big deal than it might seem! 

-   If y is actually quadratic in x, then y is linear in x\^2! 
    -   That's feature engineering.

## 2.  Constant variance

-   Or *homoscedasticity*

-   **The variance of the errors do not depend on the values of the features.**

-   Don't make bigger prediction errors for some values of x than for others.

## 3.  Normality 

-   **The errors should be independent and normally distributed.**

-   A scatter plot of target variable value and residual (model error) should look like white noise.

## 4.  Lack of perfect multicollinearity 

-   **None predictors should be a perfect linear combination of others.**

-   This can happen if you over-engineer features
    -   This is uncommon. 
    -   You'll see an error that your coefficient matrix is singular or something.

## 5.  Exogeneity 

-   **Model errors should be independent of the values of the features.**

-   In particular, errors should have mean zero. 

-   It's always good to look at a histogram of your residuals (see also normality).

## First Test

-   Determine whether the errors are normally distributed, like Shapiro-Wilk (also, plot them).

![](images/residuals.png)

## 5 Assumptions of Linear Regression: testing

-   Second I would always look at fitted value vs. residual to check homoscedasticity.

![](images/heterosceda.png)

![](images/homosceda.png)

-   For more, see for example https://people.duke.edu/\~rnau/testing.htm