---
title: "More KNN & trees"
subtitle: "Stat 253"
author: "Your Name"
output: 
  html_document:
    toc: true
    toc_float: true
---


```{r include = FALSE}
# OPTIONAL: Set a more color blind friendly palette 
palette("Okabe-Ito")
scale_colour_discrete <- function(...) scale_colour_manual(values = palette())
scale_fill_discrete   <- function(...) scale_fill_manual(values = palette())

knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env='figure',
  fig.pos = 'h',
  fig.align = 'center')
```





# Small Group Discussion

**GOAL**

Build and evaluate **nonparametric** classification models of some **categorical** outcome y.




**LAND DATA**

We'll continue to classify land type using metrics from aerial photographs:

```{r}
# Load packages
library(tidyverse)
library(tidymodels)
library(kknn)        # for KNN
library(rpart.plot)  # for plotting trees
library(class)       # for my plotting functions
```

```{r}
# Load & process data
# NOTE: Our y variable has to be converted to a "factor" variable, not a character string
land <- read.csv("https://bcheggeseth.github.io/253_spring_2024/data/land_cover.csv") %>% 
  rename(type = class) %>% 
  mutate(type = as.factor(type))
```

```{r}
# There are 9 land types! 
# Let's consider all of them (not just asphalt, grass, trees)
land %>% 
  count(type)
```

```{r}
# There are 675 data points and 147 potential predictors of land type!
dim(land)
```

```{r}
# For now: we'll consider Mean_G & NDVI
# (It's not easy to visually distinguish between 9 colors!)

# Plot type vs Mean_G alone
ggplot(land, aes(x = Mean_G, fill = type)) + 
  geom_density(alpha = 0.5) + 
  scale_fill_manual(values = c("#000000", "darkgray", "red", "#7570B3", "lightgreen", "blue", "#E6AB02", "brown", "#66A61E")) + 
  theme_minimal()
```

```{r}
# Plot type vs NDVI alone
ggplot(land, aes(x = NDVI, fill = type)) + 
  geom_density(alpha = 0.5) + 
  scale_fill_manual(values = c("#000000", "darkgray", "red", "#7570B3", "lightgreen", "blue", "#E6AB02", "brown", "#66A61E")) + 
  theme_minimal()
```

```{r}
# Finally consider both
ggplot(land, aes(y = Mean_G, x = NDVI, color = type)) + 
  geom_point() + 
  scale_color_manual(values = c("#000000", "darkgray", "red", "#7570B3", "lightgreen", "blue", "#E6AB02", "brown", "#66A61E")) + 
  theme_minimal()
```







\
\
\
\
\
\






**EXAMPLE 1: KNN**

Check out the classification regions for two KNN models of land `type` by `Mean_G` and `NDVI`: using K = 1 neighbor and using K = 16 neighbors.
Though the KNN models were *built* using standardized predictors, the predictors are plotted on their original scales here.

What questions do you have about...

- the impact of K
- how classifications are made within each neighborhood
- anything else KNN related


![](https://ajohns24.github.io/images/stat253/knn_land.png)










\
\
\
\
\
\









**EXAMPLE 2: Pruned tree**

Next, consider a **PRUNED classification tree** of land `type` by `Mean_G` and `NDVI` that was pruned as follows:

- set maximum depth to 30
- set minimum number of data points per node to 2
- tune the **cost complexity parameter**.

What questions do you have about…

- pruning trees
- using trees to make classifications
- anything else tree related

![](https://ajohns24.github.io/images/stat253/trees_land_1.png)

![](https://ajohns24.github.io/images/stat253/trees_land_2.png)







\
\
\
\
\
\





**EXAMPLE 3: Unpruned tree**

Finally, consider a (mostly) **UNPRUNED classification tree** of land `type` by `Mean_G` and `NDVI` that was built using the following tuning parameters:

- set maximum depth to 30
- set minimum number of data points per node to 2
- set cost complexity parameter to 0.


Check out the classification regions defined by this tree:

![](https://ajohns24.github.io/images/stat253/trees_land_4.png)


And the tree itself.
This tree was plotted using a function that draws the length of each branch split to be proportional to its improvement to the classification accuracy.
The labels are left off to just focus on structure:

![](https://ajohns24.github.io/images/stat253/trees_land_3.png)


- What happens to the length of the split branches the further down the tree we get? What does this mean?
- What are your thoughts about this tree?
- What questions do you have about the impact of the tuning parameters, or anything else tree related?


















 




\
\
\
\
\
\








# Exercises: Part 1

Let's focus on trees.
We've explored some general ML themes already:

- model building: parametric vs nonparametric       
    benefits and drawbacks of using nonparametric trees vs parametric algorithm logistic regression
    
- algorithm details:
    - steps of the tree algorithm
    - **tuning** a tree algorithm

In the exercises, you'll explore some other important themes:

- model building: variable selection        
    We have *147* potential predictors of land type. How can we choose which ones to use?
    
- model evaluation & comparison       
    How good are our trees? Which one should we pick?
    
- regression vs classification        
    KNN works in both settings. So do trees!



\
\



## Exercise 1: Build 10 trees

In modeling land `type` by **all 147** possible predictors, our goal will be to build a tree that optimize the `cost_complexity` parameter.

The chunk below builds 10 trees, each using a different `cost_complexity` value while fixing our other tuning parameters to be very liberal / not restrictive:        
- `tree_depth` = 30 (the default)
- `min_n` in each leaf node = 2 (the default is 20)

Reflect on the code and pause to answer any questions (Q). 

```{r}
# STEP 1: tree specification
# Q: WHAT IS NEW HERE?!
tree_spec <- decision_tree() %>%
  set_mode("classification") %>% 
  set_engine(engine = "rpart") %>% 
  set_args(cost_complexity = tune(),  
           min_n = 2, 
           tree_depth = NULL)
```

```{r}
# STEP 2: variable recipe
# NOTHING IS NEW HERE & THERE ARE NO PREPROCESSING STEPS!
variable_recipe_big <- recipe(type ~ ., data = land)

# STEP 3: tree workflow
# NOTHING IS NEW HERE!
tree_workflow_big <- workflow() %>% 
  add_recipe(variable_recipe_big) %>% 
  add_model(tree_spec)
```


```{r}
# STEP 4: Estimate 10 trees using a range of possible cost complexity values
# cost_complexity is on the log10 scale (10^(-5) to 10^(0.1))
# Q: BY WHAT CV METRIC ARE WE COMPARING THE TREES?
# Q: WHY NOT USE CV MAE?
set.seed(253)
tree_models_big <- tree_workflow_big %>% 
  tune_grid(
    grid = grid_regular(cost_complexity(range = c(-5, 0.1)), levels = 10),
    resamples = vfold_cv(land, v = 10),
    metrics = metric_set(accuracy)
  )
```






\
\



## Exercise 2: Whew!

We only built 10 trees above, and it took quite a bit of time. Why are trees computationally expensive?






\
\


## Exercise 3: Compare and finalize the tree

Just as with our other algorithms with tuning parameters, we can use the CV metrics to compare the 10 trees, and pick which one we prefer.
Here, we'll pick the *parsimonious tree* (which also happens to be the tree with the largest CV accuracy!).

Run & reflect upon the code below.

You don't have to write anything out (this is largely review), but convince yourself that you could:

- interpret the plot
- i.d. where our parsimonious tree falls on this plot
- explain what "parsimonious" means

```{r}
# Compare the CV metrics for the 10 trees
tree_models_big %>% 
  autoplot() + 
  scale_x_continuous()

# Pick the parsimonious parameter
parsimonious_param_big <- tree_models_big %>% 
  select_by_one_std_err(metric = "accuracy", desc(cost_complexity))
parsimonious_param_big

# Finalize the tree with parsimonious cost complexity
big_tree <- tree_workflow_big %>% 
  finalize_workflow(parameters = parsimonious_param_big) %>% 
  fit(data = land)
```




\
\





## Exercise 4: Examine the tree
    
Let's examine the final, tuned `big_tree` which models land `type` by all 147 predictors:       
    
```{r fig.width = 8, fig.height = 6}
# This tree has a bunch of info
# BUT: (1) it's tough to read; and (2) all branch lengths are the same (not proportional to their improvement)
big_tree %>% 
  extract_fit_engine() %>% 
  rpart.plot()
```

```{r fig.width = 8, fig.height = 6}
# This tree (1) is easier to read; and (2) plots branch lengths proportional to their improvement
# But it has less info about classification accuracy
big_tree %>% 
  extract_fit_engine() %>%
  plot()
big_tree %>% 
  extract_fit_engine() %>%
  text(cex = 0.8)
```

Use the *second* tree plot to answer the following questions.

a. Are any land types *not* captured somewhere in the tree: asphalt, building, car, concrete, grass, pool, shadow, soil, tree? If so, why do you think this is?




b. This tree considered all 147 possible predictors. Do all of these appear in the final tree?




c. Of the predictors used in the tree, which seem to be the most important?







    
    
\
\



## Exercise 5: Identifying useful predictors
    
Luckily, we don't have to *guesstimate* the importance of different predictors.

There's a mathematical metric: **variable importance**.

Roughly, a predictor's importance measures the total improvement in node purity if we were to split on this predictor (even if the predictor isn't ultimately used in a split).

The bigger the better!

Check out the importance of our 147 predictors.

```{r}
# Check out just the 10 most important (for simplicity)
big_tree %>%
  extract_fit_engine() %>%
  pluck("variable.importance") %>% 
  head(10)
```

And plot these to help with comparison:

```{r fig.width = 8, fig.height = 8}
# By default, this plots only 10 predictors
# At most, it will plot only half of our predictors here
# (I'm not sure what the max is in general!)
library(vip)
big_tree %>% 
  vip(geom = "point", num_features = 147)
```

a. What are the 3 most important predictors by this metric?




b. Do these 3 predictors appear in our tree?




c. *Why* do you think a predictor can have high importance but not appear in the tree? Name 2 reasons.









\
\



## Exercise 6: How good is the tree?!? Part 1

Stepping back, our goal for building this tree was to classify land type in pixels in images.
As an example, suppose we have an image like the first one in our data set:

```{r}
head(land, 1)
```

Our tree correctly predicts that this is a car:

```{r}
big_tree %>% 
  predict(new_data = head(land, 1))
```

But how good are the classifications overall?
Let's first consider the CV *overall accuracy rate*:   

```{r}
parsimonious_param_big
```

a. *Interpret* 0.804, the CV overall accuracy rate for the `big_tree`.




b. Calculate & compare the tree's CV overall accuracy rate to the **no information rate**: is the `big_tree` better than just always guessing the most common land type? You'll need the following info:     

```{r}
land %>% 
  nrow()
  
land %>% 
  count(type)
```



c. Why can't we calculate, thus compare, the sensitivity & specificity of our tree? HINT: Think of how these are defined.
    
    
    
    
    
\
\




## Exercise 7: How good is the tree?!? Part 2

The above CV metric gives us a sense of the *overall* quality of using our tree to classify *new* pixels in images.

But it doesn't give any insight into the quality of classifications for any particular land type.

To that end, let's consider the *in-sample* confusion matrix (i.e. how well our tree classified the pixels in images in our sample).

NOTE: We could also get a CV version, but the code is long and the CV accuracy will do!

```{r fig.height = 6, fig.width = 8}
in_sample_confusion <- big_tree %>% 
  augment(new_data = land) %>% 
  conf_mat(truth = type, estimate = .pred_class)
in_sample_confusion

# The mosaic plot of this matrix is too messy to be very useful here
in_sample_confusion %>% 
  autoplot() +    
  aes(fill = rep(colnames(in_sample_confusion$table), ncol(in_sample_confusion$table))) + 
  scale_fill_manual(values = c("#000000", "darkgray", "red", "#7570B3", "lightgreen", "blue", "#E6AB02", "brown", "#66A61E")) + 
  theme(legend.position = "none")
```    

a. Confirm that 96.6% of asphalt pixels in images were correctly classified as asphalt.

```{r}

```

b. What land type was the *hardest* for the tree to classify? Why might this be?







\
\
\
\
\
\


# Exercises: Part 2

Just like KNN, trees can be applied in both regression *and* classification settings.
Thus trees add to our collection of *nonparametric* regression techniques, including KNN, LOESS, and GAM.

To explore, we'll use regression trees to model the `body_mass_g` of penguins by their `bill_length_mm` and `species`.

*This is for demonstration purposes only!!!*
As the plot below demonstrates, this relationship isn't complicated enough to justify using a nonparametric algorithm:

```{r}
data(penguins)
ggplot(penguins, aes(y = body_mass_g, x = bill_length_mm, color = species)) + 
  geom_point()
```


Run the following code to build a regression tree of this relationship.
Pause to reflect upon the questions in the comments:

```{r}
# CHUNK GOAL
# Build a bunch of trees using different cost complexity parameters

# STEP 1: regression tree specification
# QUESTION: How does this differ from our classification tree specification?
tree_spec <- decision_tree() %>%
  set_mode("regression") %>% 
  set_engine(engine = "rpart") %>% 
  set_args(cost_complexity = tune(),  
           min_n = 2, 
           tree_depth = 20)

# STEP 2: variable recipe
# NOTHING IS NEW HERE!
variable_recipe <- recipe(body_mass_g ~ bill_length_mm + species, data = penguins)

# STEP 3: tree workflow
# NOTHING IS NEW HERE!
tree_workflow <- workflow() %>% 
  add_recipe(variable_recipe) %>% 
  add_model(tree_spec)

# STEP 4: Estimate multiple trees using a range of possible cost complexity values
# QUESTION: How do the CV metrics differ from our classification tree?
set.seed(253)
tree_models <- tree_workflow %>% 
  tune_grid(
    grid = grid_regular(cost_complexity(range = c(-5, -1)), levels = 10),
    resamples = vfold_cv(penguins, v = 10),
    metrics = metric_set(mae)
  )
```

```{r}
# CHUNK GOAL:
# Finalize the tree using the parsimonious cost complexity parameter
# NOTHING IS NEW HERE

# Identify the parsimonious cost complexityparameter
parsimonious_param <- tree_models %>% 
  select_by_one_std_err(metric = "mae", desc(cost_complexity))

# Finalize the tree with parsimonious cost complexity
regression_tree <- tree_workflow %>% 
  finalize_workflow(parameters = parsimonious_param) %>% 
  fit(data = penguins)
```




\
\


## Exercise 8: Regression tree

Check out the resulting regression tree:

```{r}
# This code is the same as for classification trees!
regression_tree %>% 
  extract_fit_engine() %>% 
  rpart.plot()
```

a. Use your tree (by "hand") to predict the body mass for the 2 following penguins:       

```{r}
new_penguins <- data.frame(species = c("Adelie", "Gentoo"), bill_length_mm = c(45, 45))
new_penguins
```



b. Check your work in part a:       

```{r}
regression_tree %>% 
  augment(new_data = new_penguins)
```



c. Regression trees partition the data points into separate prediction regions. Check out this tree's predictions (dark dots) and convince yourself that these are consistent with the tree:    

```{r fig.width = 6}
regression_tree %>% 
  augment(new_data = penguins) %>% 
  ggplot(aes(x = bill_length_mm, y = body_mass_g, color = species)) + 
    geom_point(alpha = 0.35, size = 0.5) + 
    geom_point(aes(x = bill_length_mm, y = .pred, color = species), size = 1.5) +
    facet_wrap(~ species) + 
    theme_minimal()
```
  
d. Based on what you've observed here, what do you think is a drawback of regression trees?   






\
\




## Exercise 9: Visual essay

Check out this [visual essay](http://www.r2d3.us/visual-intro-to-machine-learning-part-1/) on trees!



\
\



## Exercise 10: Work on Homework

With your remaining time, work on HW 5 (logistic regression) or HW 6 (KNN, trees)
    
    
   



# R code: KNN

Suppose we want to build a KNN model of some categorical response variable y using predictors x1 and x2 in our sample_data.

```{r eval = FALSE}
# Load packages
library(tidymodels)
library(kknn)

# Resolves package conflicts by preferring tidymodels functions
tidymodels_prefer()
```


\
\


**Build models for a variety of tuning parameters K**


```{r eval = FALSE}
# STEP 1: KNN model specification
knn_spec <- nearest_neighbor() %>%
  set_mode("classification") %>% 
  set_engine(engine = "kknn") %>% 
  set_args(neighbors = tune())

# STEP 2: variable recipe
variable_recipe <- recipe(y ~ x1 + x2, data = sample_data) %>% 
  step_nzv(all_predictors()) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  step_normalize(all_numeric_predictors())

# STEP 3: KNN workflow
knn_workflow <- workflow() %>% 
  add_recipe(variable_recipe) %>% 
  add_model(knn_spec)

# STEP 4: Estimate multiple KNN models using a range of possible K values
set.seed(___)
knn_models <- knn_workflow %>% 
  tune_grid(
    grid = grid_regular(neighbors(range = c(___, ___)), levels = ___),
    resamples = vfold_cv(sample_data, v = ___),
    metrics = metric_set(accuracy)
  )
```



\
\


**Tuning K**

```{r eval = FALSE}
# Calculate CV accuracy for each KNN model
knn_models %>% 
  collect_metrics()

# Plot CV accuracy (y-axis) for the KNN model from each K (x-axis)
knn_models %>% 
  autoplot()

# Identify K which produced the highest (best) CV accuracy
best_K <- select_best(knn_models, metric = "accuracy")
best_K

# Get the CV accuracy for KNN when using best_K
knn_models %>% 
  collect_metrics() %>% 
  filter(neighbors == best_K$neighbors)
```


\
\


**Finalizing the "best" KNN model**

```{r eval = FALSE}
final_knn_model <- knn_workflow %>% 
  finalize_workflow(parameters = best_k) %>% 
  fit(data = sample_data)
```


\
\


**Use the KNN to make predictions / classifications**


```{r eval = FALSE}
# Put in a data.frame object with x1 and x2 values (at minimum)
final_knn_model %>% 
  predict(new_data = ___)  
```






\
\
\
\
\
\



# R code: trees

Suppose we want to build a tree of some categorical response variable y using predictors x1 and x2 in our sample_data.

```{r eval = FALSE}
# Load packages
library(tidymodels)
library(rpart)
library(rpart.plot)

# Resolves package conflicts by preferring tidymodels functions
tidymodels_prefer()
```


\
\


**Build trees for a variety of tuning parameters**

We'll focus on optimizing the `cost_complexity` parameter, while setting `min_n` and `tree_depth` to fixed numbers.
For example, setting `min_n` to 2 and `tree_depth` to 30 set only loose restrictions, letting `cost_complexity` do the pruning work.



```{r eval = FALSE}
# STEP 1: tree specification
# If y is quantitative, change "classification" to "regression"
tree_spec <- decision_tree() %>%
  set_mode("classification") %>% 
  set_engine(engine = "rpart") %>% 
  set_args(cost_complexity = tune(),  
           min_n = 2, 
           tree_depth = 30)

# STEP 2: variable recipe
# There are no necessary preprocessing steps for trees!
variable_recipe <- recipe(y ~ x1 + x2, data = sample_data)

# STEP 3: tree workflow
tree_workflow <- workflow() %>% 
  add_recipe(variable_recipe) %>% 
  add_model(tree_spec)

# STEP 4: Estimate multiple trees using a range of possible cost complexity values
# - If y is quantitative, change "accuracy" to "mae"
# - cost_complexity is on the log10 scale (10^(-5) to 10^(0.1))
#   I start with a range from -5 to 2 and then tweak
set.seed(___)
tree_models <- tree_workflow %>% 
  tune_grid(
    grid = grid_regular(cost_complexity(range = c(___, ___)), levels = ___),
    resamples = vfold_cv(sample_data, v = ___),
    metrics = metric_set(accuracy)
  )
```

\
\


**Tuning cost complexity**

```{r eval = FALSE}
# Plot the CV accuracy vs cost complexity for our trees
# x-axis is on the original (not log10) scale
tree_models %>% 
  autoplot() + 
  scale_x_continuous()

# Identify cost complexity which produced the highest CV accuracy
best_cost <- tree_models %>% 
  select_best(metric = "accuracy")

# Get the CV accuracy when using best_cost
tree_models %>% 
  collect_metrics() %>% 
  filter(cost_complexity == best_cost$cost_complexity)

# Identify cost complexity which produced the parsimonious tree
parsimonious_cost <- tree_models %>% 
  select_by_one_std_err(metric = "accuracy", desc(cost_complexity))
```



\
\



**Finalizing the tree**

```{r eval = FALSE}
# Plug in best_cost or parsimonious_cost
final_tree <- tree_workflow %>% 
  finalize_workflow(parameters = ___) %>% 
  fit(data = sample_data)
```


\
\


**Plot the tree**

```{r eval = FALSE}
# Tree with accuracy info in each node
# Branches are NOT proportional to classification improvement
final_tree %>% 
  extract_fit_engine() %>% 
  rpart.plot()

# Tree withOUT accuracy info in each node
# Branches ARE proportional to classification improvement
final_tree %>% 
  extract_fit_engine() %>% 
  plot()
final_tree %>% 
  extract_fit_engine() %>% 
  text()
```




\
\


**Use the tree to make predictions / classifications**


```{r eval = FALSE}
# Put in a data.frame object with x1 and x2 values (at minimum)
final_tree %>% 
  predict(new_data = ___)  

# OR 
final_tree %>% 
  augment(new_data = ___)
```




\
\


**Examine variable importance**

```{r eval = FALSE}
# Print the metrics
final_tree %>%
  extract_fit_engine() %>%
  pluck("variable.importance")

# Plot the metrics
# Plug in the number of top predictors you wish to plot
# (The upper limit varies by application!)
library(vip)
final_tree %>% 
  vip(geom = "point", num_features = ___)
```



\
\



**Evaluate the classifications using in-sample metrics**

```{r eval = FALSE}
# Get the in-sample confusion matrix
in_sample_confusion <- final_tree %>% 
  augment(new_data = sample_data) %>% 
  conf_mat(truth = type, estimate = .pred_class)
in_sample_confusion

# Plot the matrix using a mosaic plot
# See exercise for what to do when there are more categories than colors in our pallette!
in_sample_confusion %>% 
  autoplot() +    
  aes(fill = rep(colnames(in_sample_confusion$table), ncol(in_sample_confusion$table))) + 
  theme(legend.position = "none")
```




\
\
\
\
\
\


