library(tidyverse)
library(class)
library(ISLR)
library(plotly) # 3D Graphs
library(MASS) # LDA

wine_full <- read.csv("/Users/Will/Documents/College/Junior/Spring 2024/STAT218/Datasets/Wine Training.csv") 
wine_test <- read.csv("/Users/Will/Documents/College/Junior/Spring 2024/STAT218/Datasets/Will_wine_data.csv")

To determine which 3 features we want to select and use to make predictions, we will graph the “quality” compared to all of the features to look for relationships. A multiple regression can also point us in the right direction on correlations between features and quality.

# relationship exploration

wine_full$quality = as.character(wine_full$quality)
class(wine_full$quality)
## [1] "character"
model <- lm(quality ~ fixed.acidity + volatile.acidity + citric.acid + residual.sugar + chlorides + free.sulfur.dioxide + total.sulfur.dioxide + density + pH + sulphates + alcohol, data = wine_full)

summary(model)
## 
## Call:
## lm(formula = quality ~ fixed.acidity + volatile.acidity + citric.acid + 
##     residual.sugar + chlorides + free.sulfur.dioxide + total.sulfur.dioxide + 
##     density + pH + sulphates + alcohol, data = wine_full)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.70735 -0.37343 -0.05225  0.44159  1.93288 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.754e+01  2.535e+01   1.481   0.1389    
## fixed.acidity         3.141e-02  3.124e-02   1.005   0.3150    
## volatile.acidity     -1.218e+00  1.437e-01  -8.478  < 2e-16 ***
## citric.acid          -1.934e-01  1.783e-01  -1.085   0.2781    
## residual.sugar        1.949e-02  1.843e-02   1.058   0.2903    
## chlorides            -2.153e+00  5.453e-01  -3.949 8.36e-05 ***
## free.sulfur.dioxide   4.404e-03  2.607e-03   1.689   0.0915 .  
## total.sulfur.dioxide -3.433e-03  8.767e-04  -3.916 9.57e-05 ***
## density              -3.335e+01  2.588e+01  -1.288   0.1979    
## pH                   -3.023e-01  2.292e-01  -1.319   0.1875    
## sulphates             9.520e-01  1.329e-01   7.164 1.44e-12 ***
## alcohol               2.277e-01  3.108e-02   7.328 4.55e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6373 on 1087 degrees of freedom
## Multiple R-squared:  0.3619, Adjusted R-squared:  0.3554 
## F-statistic: 56.04 on 11 and 1087 DF,  p-value: < 2.2e-16
wine_full |>
  plot_ly(x = ~sulphates, y = ~volatile.acidity, z = ~alcohol, color = ~quality, colors = c("red", "red", "darkgrey", "darkgrey", "steelblue", "steelblue")) |> 
  add_markers()

From these diagrams, we see that there appear to be strong relationships between quality and volatile acidity, alcohol, and sulphates.

Algorithm Selection:

## 
## high  low midd 
##  143   38  918

Given that our data is distributed rather unevenly, kNN here might be misleading due to the high number of quality 5’s and 6’s. We want to find the extremes. Therefore, an LDA model would likely perform better as it will take into account the size of these groups.

# splitting our data into test and train
rows_to_keep <- sample(1:nrow(wine_full), size = nrow(wine_full)/2)
wine_full_test <- wine_full[rows_to_keep,]
wine_full_train <- wine_full[-rows_to_keep,]
# LDA
lda <- lda(quality ~ sulphates + alcohol + volatile.acidity, data = wine_full_train)

preds <- predict(lda, wine_full_test)

table(wine_full_test$quality, preds$class)
##    
##       3   4   5   6   7   8
##   3   1   0   3   0   0   0
##   4   0   0  14   3   1   0
##   5   2   0 171  49   3   0
##   6   0   0  94 115  19   0
##   7   0   0   5  41  24   0
##   8   0   0   0   2   2   0

To evaluate this model, we really only care about how well it performs at the extremes, that is low quality groups and high quality groups. And it’s not that good at it given our confusion matrix above… Let’s try to change our model to only predict middle, low, and high quality instead of exact scores.

# LDA 2 Groups
lda_groups <- lda(quality_group ~ sulphates + alcohol + volatile.acidity, data = wine_full_train)

pred_probs_groups <- predict(lda_groups, wine_full_test)

table(wine_full_test$quality_group, pred_probs_groups$class)
##       
##        high low midd
##   high   25   0   49
##   low     1   1   20
##   midd   21   1  431
# accuracy measurement from straight up LDA class prediction
eval <- data.frame(pred_probs_groups, wine_full_test$quality_group)

bottomAccuracy <- eval[order(eval$posterior.low, decreasing = TRUE), ][1:10,] |>
  mutate(success = if_else(class == "low" & wine_full_test.quality_group == "low", 1, 0)) |>
  pull(success) |>
  sum() / 10
           
bottomAccuracy
## [1] 0.1
topAccuracy <- eval[order(eval$posterior.high, decreasing = TRUE), ][1:10,] |>
  mutate(success = if_else(class == "high" & wine_full_test.quality_group == "high", 1, 0)) |>
  pull(success) |>
  sum() / 10

topAccuracy
## [1] 0.5
# LDA Tweaking with Accuracies

lda_data <- data.frame(preds)

lda_data <- lda_data |> 
  mutate(high = 10*posterior.8 + posterior.7 - posterior.6) |> 
  mutate(low = 10*posterior.3 + posterior.4 - posterior.5) |> 
  cbind(wine_full_test)

top_accuracy <- lda_data[order(lda_data$high, decreasing = TRUE), ][1:10, ] |> 
  mutate(success = if_else(quality_group == "high", 1, 0)) |> 
  pull(success) |> 
  sum() / 10

top_accuracy
## [1] 0.5
bottom_accuracy <- lda_data[order(lda_data$low, decreasing = TRUE), ][1:10, ] |> 
  mutate(success = if_else(quality_group == "low", 1, 0)) |> 
  pull(success) |> 
  sum() / 10

bottom_accuracy
## [1] 0.4