Coffee Quality Score Modelling

Following on from the exploratory data analysis of reviews of coffee quality from different suppliers, we will now attempt to model the main quality metric, Total.Cup.Points, based on the characteristics of the beans used to make the coffee.

We have 1338 rows of data, of which the first 10 rows are given below.

We would like to develop a random forest model to describe how Total.Cup.Points depends on the input variables of this data. To choose the optimal parameters to use with our model, we would like to crossvalidate using training data, before testing the chosen parameters using test data.

Feature Importance

To estimate the relative importance of the features, we fit an unoptimised random forest model to all the data and determine variable importance from the model.

The feature variables with greatest influence on quality appear to be:

• Altitude
• Category.Two.Defects
• Owner
• Moisture
• Weight of Sample Bag
• Grading Month
• Variety
• Category.One.Defects
• Harvest and Grading Year

Processing method, Quakers, Species and Color appear to not be as important.

We have difficulties with this data set in that some of the categorical variables have low frequency values within them. These factor levels may not appear in all folds of the crossvalidation data, and may not be present in both the training and test data. If we use our workflow to just create dummies for these variables as is, we will encounter errors during prediction steps when the training and testing variables don’t match.

To compensate for this, we will identify our top 14 countries of origin, top 6 Owners, top 5 varieties and the top 3 processing methods and change the other values in these categories to “Other”. We are reducing the number of distinct categories for each of these variables.

After removing rows with missing data for Owner, Country or Variety, a random 20% will be removed to act as the test set. We will also randomly define the 4-fold crossvalidation splits for the training data.

Training data rows = 890 
Testing data rows = 222 

Set up Pre-treatment Recipes

We now define the pretreatment recipes. We are going to impute some data into missing values for crossvalidation. We will impute the mean altitude, the median Harvest.Year, and the mode for the categorical variables processing method and Color. We will also normalize numeric variables and remove those that are highly correlated. We will create variables for grading month and year and create dummy variables for the remaining categorical variables.

recipe_1 <-
  train %>%
  recipe(Total.Cup.Points ~ .) %>%
  step_meanimpute(mean_alt) %>%
  step_medianimpute(Harvest.Year) %>%
  step_modeimpute(pm, Color) %>%
  step_normalize(all_numeric(), -all_outcomes(), -Harvest.Year) %>%
  step_date(new_date, features = c("month", "year")) %>%
  step_rm(new_date) %>% 
  step_dummy(Owner, Species, Color, coo, var, pm, new_date_month) %>% 
  step_corr(all_numeric(), -all_outcomes())

What does the subsequent training frame look like?

Random Forest Model

We will create a tuning workflow in order to optimise hyperparameters for a random forest regression model.

We will be optimising trees and mtry values for our ranger model, using crossvalidation results of our training data. Mean absolute error (MAE) will be evaluated for each set of predictions.

rf_model <- 
  rand_forest(trees = tune(), mtry = tune()) %>% 
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("regression")

grid_rf <- 
  grid_max_entropy(
    trees(range = c(100, 2000)),
    mtry(range = c(3, 15)),
    size = 30)

my_metrics <- metric_set(mae)

After crossvalidating for a range of trees and mtry values, let’s plot values of MAE and list the best models.

We will pick the values of trees and mtry with the lowest error. Now we will choose these values to train a model using the training data and fit the model to the test data. Let’s see how the predictions compare with reality for the test set and calculate the test MAE. The line on the plot is the line of equality.

Mean Absolute Error for Test set = 1.421002

Our test set error is better than our best crossvalidation error.

Linear Model

Let’s compare this with a linear regression model. We will be optimising mixture and penalty values for our glmnet model, using crossvalidation results of our training data. Mean absolute error (MAE) will be evaluated for each set of predictions.

model_glm <- 
  linear_reg(penalty = tune(), mixture = tune()) %>% 
  set_engine("glmnet")

grid_glm <- 
  grid_max_entropy(
    penalty(range = c(0.01, 0.3), trans = NULL),
    mixture(range = c(0, 1)),
    size = 30)

my_metrics <- metric_set(mae)

After crossvalidating for a range of mixture and penalty values, let’s plot values of MAE and list the best models.

Again we will pick the values of mixture and penalty with the lowest error. Now we will choose these values to train a model using the training data and fit the model to the test data. Again let’s see how the predictions compare with reality for the test set and calculate the test MAE. The line on the plot is the line of equality.

Mean Absolute Error for Test set = 1.479737

Again, the error on the test set is better than that of the crossvalidation results.

Overall, it appears that the random forest model is better at predicting coffee quality than the linear regression model. It may be possible to reduce the prediction error further by combining the two model predictions.

---

End