knitr::opts_chunk$set(message = FALSE, warning = FALSE)
options(cli.unicode = FALSE)
library(mlbench)
library(tidyverse)
library(e1071)
library(caret)
library(corrplot)

Exercise 3.1

The UC Irvine Machine Learning Repository contains a data set related to glass identification. The data consist of 214 glass samples labeled as one of seven class categories. There are nine predictors: the refractive index and the percentages of eight elements (Na, Mg, Al, Si, K, Ca, Ba and Fe).

data(Glass)
str(Glass)
## 'data.frame':    214 obs. of  10 variables:
##  $ RI  : num  1.52 1.52 1.52 1.52 1.52 ...
##  $ Na  : num  13.6 13.9 13.5 13.2 13.3 ...
##  $ Mg  : num  4.49 3.6 3.55 3.69 3.62 3.61 3.6 3.61 3.58 3.6 ...
##  $ Al  : num  1.1 1.36 1.54 1.29 1.24 1.62 1.14 1.05 1.37 1.36 ...
##  $ Si  : num  71.8 72.7 73 72.6 73.1 ...
##  $ K   : num  0.06 0.48 0.39 0.57 0.55 0.64 0.58 0.57 0.56 0.57 ...
##  $ Ca  : num  8.75 7.83 7.78 8.22 8.07 8.07 8.17 8.24 8.3 8.4 ...
##  $ Ba  : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Fe  : num  0 0 0 0 0 0.26 0 0 0 0.11 ...
##  $ Type: Factor w/ 6 levels "1","2","3","5",..: 1 1 1 1 1 1 1 1 1 1 ...

a. Using visualizations, explore the predictor variables to understand their distributions as well as the relationships between predictors.

glass_long <- Glass |>
  select(-Type) |>
  pivot_longer(everything(), names_to = "predictor", values_to = "value")

ggplot(glass_long, aes(value)) +
  geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
  facet_wrap(~ predictor, scales = "free") +
  labs(title = "Distributions of the nine predictors")

ggplot(glass_long, aes(y = value)) +
  geom_boxplot(fill = "steelblue") +
  facet_wrap(~ predictor, scales = "free") +
  labs(title = "Boxplots of the nine predictors") +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

summary(Glass[, 1:9])
##        RI              Na              Mg              Al       
##  Min.   :1.511   Min.   :10.73   Min.   :0.000   Min.   :0.290  
##  1st Qu.:1.517   1st Qu.:12.91   1st Qu.:2.115   1st Qu.:1.190  
##  Median :1.518   Median :13.30   Median :3.480   Median :1.360  
##  Mean   :1.518   Mean   :13.41   Mean   :2.685   Mean   :1.445  
##  3rd Qu.:1.519   3rd Qu.:13.82   3rd Qu.:3.600   3rd Qu.:1.630  
##  Max.   :1.534   Max.   :17.38   Max.   :4.490   Max.   :3.500  
##        Si              K                Ca               Ba       
##  Min.   :69.81   Min.   :0.0000   Min.   : 5.430   Min.   :0.000  
##  1st Qu.:72.28   1st Qu.:0.1225   1st Qu.: 8.240   1st Qu.:0.000  
##  Median :72.79   Median :0.5550   Median : 8.600   Median :0.000  
##  Mean   :72.65   Mean   :0.4971   Mean   : 8.957   Mean   :0.175  
##  3rd Qu.:73.09   3rd Qu.:0.6100   3rd Qu.: 9.172   3rd Qu.:0.000  
##  Max.   :75.41   Max.   :6.2100   Max.   :16.190   Max.   :3.150  
##        Fe         
##  Min.   :0.00000  
##  1st Qu.:0.00000  
##  Median :0.00000  
##  Mean   :0.05701  
##  3rd Qu.:0.10000  
##  Max.   :0.51000

Relationships between the predictors, using the same kind of correlation matrix the chapter uses for the cell segmentation data:

corrplot(cor(Glass[, 1:9]), method = "number", type = "upper",
         tl.col = "black", number.cex = 0.7)

The distributions are very different from one another:

  • Na, Si and Al are roughly symmetric and bell shaped.
  • RI, Ca, K and Fe are right-skewed with long upper tails.
  • Mg is left-skewed and bimodal, with a large group of samples at 0 and another group near 3.5.
  • Ba and Fe are zero in most samples, which is the near-zero variance pattern described in Sect. 3.5.

Most pairs of predictors are only weakly correlated, so collinearity is much less of a problem here than in the chapter’s cell segmentation example. The important exception is RI and Ca, with a correlation of 0.81, which makes sense physically, since calcium content drives the refractive index. There are also moderate relationships between RI and Si (-0.54), Mg and Al (-0.48), Mg and Ba (-0.49), and Al and Ba (0.48).

table(Glass$Type)
## 
##  1  2  3  5  6  7 
## 70 76 17 13  9 29

The exercise describes seven class categories, but the data set in mlbench has six: there are no samples of type 4. The classes are also unbalanced, since types 1 and 2 hold two thirds of the samples while type 6 has only 9.

b. Do there appear to be any outliers in the data? Are any predictors skewed?

The chapter gives two diagnostics for skewness: the ratio of the largest value to the smallest value (a ratio above about 20 indicates significant skewness), and the skewness statistic, which is near zero for a symmetric distribution, positive for right skew and negative for left skew.

outlier_count <- function(x) {
  q <- quantile(x, c(0.25, 0.75))
  sum(x < q[1] - 1.5 * IQR(x) | x > q[2] + 1.5 * IQR(x))
}

tibble(
  predictor = names(Glass)[1:9],
  min = round(sapply(Glass[, 1:9], min), 2),
  max = round(sapply(Glass[, 1:9], max), 2),
  max_min_ratio = round(sapply(Glass[, 1:9],
                               function(x) ifelse(min(x) > 0, max(x) / min(x), NA)), 1),
  skewness = round(sapply(Glass[, 1:9], skewness), 2),
  outliers = sapply(Glass[, 1:9], outlier_count),
  zeros = sapply(Glass[, 1:9], function(x) sum(x == 0))
)
## # A tibble: 9 x 7
##   predictor   min   max max_min_ratio skewness outliers zeros
##   <chr>     <dbl> <dbl>         <dbl>    <dbl>    <int> <int>
## 1 RI         1.51  1.53           1       1.6        17     0
## 2 Na        10.7  17.4            1.6     0.45        7     0
## 3 Mg         0     4.49          NA      -1.14        0    42
## 4 Al         0.29  3.5           12.1     0.89       18     0
## 5 Si        69.8  75.4            1.1    -0.72       12     0
## 6 K          0     6.21          NA       6.46        7    30
## 7 Ca         5.43 16.2            3       2.02       26     0
## 8 Ba         0     3.15          NA       3.37       38   176
## 9 Fe         0     0.51          NA       1.73       12   144

Outliers. Yes. Every predictor except Mg has values outside the usual 1.5 × IQR range. Ba has the most (38) and Ca has 26. Most of the Ba and Fe flags are not really errors, though: those two elements are zero in most samples, so any sample that actually contains them gets flagged. The clearest single outlier is a potassium reading of 6.21 when the median K is only 0.56. Ca also reaches 16.19 against a median of 8.60. As the chapter cautions, these look like a valid but poorly sampled part of the population (the barium-rich headlamp glass, for example) rather than mistakes, so they should not simply be deleted.

Skewness. Yes, several predictors are skewed, though the two diagnostics disagree here. The max-to-min ratio does not reach 20 for any predictor, and cannot even be computed for Mg, K, Ba and Fe, whose minimum is zero. That rule works poorly for these predictors, since they are percentages measured around a high baseline (Si never falls below 69.8). The skewness statistic is the more useful diagnostic:

  • K (6.5), Ba (3.4), Ca (2.0), Fe (1.7) and RI (1.6) are right-skewed.
  • Mg (-1.1) and Si (-0.7) are left-skewed.
  • Na (0.4) and Al (0.9) are close to symmetric.

Skewness beyond about ±1 is worth correcting, so K, Ba, Ca, Fe, RI and Mg all qualify.

c. Are there any relevant transformations of one or more predictors that might improve the classification model?

The chapter’s Box-Cox procedure is applied only to predictors whose values are all greater than zero, and Mg, K, Ba and Fe all contain zeros. The Yeo-Johnson transformation does the same job but allows zeros, so I use it here, together with the centering and scaling the chapter recommends before any distance-based method.

trans <- preProcess(Glass[, 1:9], method = c("YeoJohnson", "center", "scale"))
trans
## Created from 214 samples and 9 variables
## 
## Pre-processing:
##   - centered (9)
##   - ignored (0)
##   - scaled (9)
##   - Yeo-Johnson transformation (5)
## 
## Lambda estimates for Yeo-Johnson transformation:
## -0.19, 2.24, 0, -0.98, -1.36
glass_trans <- predict(trans, Glass[, 1:9])

tibble(
  predictor = names(Glass)[1:9],
  skew_before = round(sapply(Glass[, 1:9], skewness), 2),
  skew_after = round(sapply(glass_trans, skewness), 2)
)
## # A tibble: 9 x 3
##   predictor skew_before skew_after
##   <chr>           <dbl>      <dbl>
## 1 RI               1.6        1.6 
## 2 Na               0.45      -0.01
## 3 Mg              -1.14      -0.88
## 4 Al               0.89       0   
## 5 Si              -0.72      -0.72
## 6 K                6.46      -0.07
## 7 Ca               2.02      -0.21
## 8 Ba               3.37       3.37
## 9 Fe               1.73       1.73
glass_trans |>
  pivot_longer(everything(), names_to = "predictor", values_to = "value") |>
  ggplot(aes(value)) +
  geom_histogram(bins = 30, fill = "darkseagreen4", colour = "white") +
  facet_wrap(~ predictor, scales = "free") +
  labs(title = "Distributions after the Yeo-Johnson transformation")

The transformations that would help are:

  1. Yeo-Johnson (Box-Cox where values allow it), plus centering and scaling. This fixes the worst skewness: RI drops from 1.6 to about 0.3, Ca from 2.0 to about 0.4, and Al and Na become nearly symmetric. Centering and scaling matter because the predictors are on completely different scales, which affects any method based on distances, such as k-nearest neighbours or support vector machines.

  2. The transformation does not help Ba, Fe, Mg or K much. Their problem is not skewness but a pile-up of values at zero, and no single transformation can spread out a spike at one value. A better option is to replace them with a binary indicator such as “contains barium, yes or no”, or to drop Ba and Fe under the Sect. 3.5 rule for near-zero variance predictors.

  3. The spatial sign transformation (Sect. 3.3), applied after centering and scaling, would pull extreme values such as the K reading of 6.21 toward the rest of the data instead of discarding those samples. It transforms the predictors as a group, so it suits this data set, where the outliers appear in several elements at once.

  4. PCA would handle the strong RI and Ca correlation by replacing the correlated predictors with uncorrelated components. Alternatively, one of the two could be removed, since they measure much the same thing.

Exercise 3.2

The soybean data can also be found at the UC Irvine Machine Learning Repository. Data were collected to predict disease in 683 soybeans. The 35 predictors are mostly categorical and include information on the environmental conditions (such as temperature and precipitation) and the plant conditions (such as leaf spots and mold growth). The outcome has 19 distinct classes.

data(Soybean)
dim(Soybean)
## [1] 683  36

a. Investigate the frequency distributions for the categorical predictors. Are any of the distributions degenerate in the ways discussed earlier in this chapter?

Soybean |>
  select(-Class) |>
  mutate(across(everything(), as.character)) |>
  pivot_longer(everything(), names_to = "predictor", values_to = "value") |>
  ggplot(aes(value)) +
  geom_bar(fill = "steelblue") +
  facet_wrap(~ predictor, scales = "free", ncol = 5) +
  labs(title = "Frequency distributions of the 35 categorical predictors",
       x = NULL, y = "Count")

Section 3.5 gives two criteria for a near-zero variance predictor: the fraction of unique values over the sample size is low (say 10%), and the ratio of the frequency of the most prevalent value to the frequency of the second most prevalent value is large (say around 20). nearZeroVar() applies both.

nzv <- nearZeroVar(Soybean[, -1], saveMetrics = TRUE)
nzv[nzv$nzv, ]
##           freqRatio percentUnique zeroVar  nzv
## leaf.mild     26.75     0.4392387   FALSE TRUE
## mycelium     106.50     0.2928258   FALSE TRUE
## sclerotia     31.25     0.2928258   FALSE TRUE

Yes, three predictors are degenerate in exactly that sense:

table(Soybean$mycelium, useNA = "ifany")
## 
##    0    1 <NA> 
##  639    6   38
table(Soybean$sclerotia, useNA = "ifany")
## 
##    0    1 <NA> 
##  625   20   38
table(Soybean$leaf.mild, useNA = "ifany")
## 
##    0    1    2 <NA> 
##  535   20   20  108
  • mycelium: 639 samples take the value 0 and only 6 take the value 1, a frequency ratio of about 107 to 1, with 0.3% unique values.
  • sclerotia: 625 zeros against 20 ones, a ratio of about 31 to 1.
  • leaf.mild: 535 zeros against 20 and 20 in the other two categories, a ratio of about 27 to 1.

All three clear the “around 20” threshold for the frequency ratio and have well under 10% unique values. They carry almost no information for most samples, and as the chapter warns, a resampled training set could easily contain only one of their values, which would leave a zero variance predictor in that fold. Several others, such as shriveling, int.discolor, lodging and leaf.malf, are nearly as unbalanced but fall just below the cutoff.

c. Develop a strategy for handling missing data, either by eliminating predictors or imputation.

# Step 1: drop the three degenerate predictors
soybean_reduced <- Soybean |>
  select(-mycelium, -sclerotia, -leaf.mild)

# Step 2: give the remaining factors an explicit "missing" category
soybean_clean <- soybean_reduced |>
  mutate(across(-Class, ~ fct_na_value_to_level(., level = "missing")))

sum(is.na(soybean_clean))
## [1] 0
table(soybean_clean$hail)
## 
##       0       1 missing 
##     435     127     121

My strategy has three parts.

1. Remove the degenerate predictors. mycelium, sclerotia and leaf.mild are dropped for the reasons given in part a. leaf.mild is also one of the predictors with the most missing values, so this removes part of the missing data problem at the same time.

2. Do not delete the incomplete samples. The chapter notes that removing samples is reasonable for large data sets when the missingness is not informative, but costly in small ones. Both warnings apply here. Deleting would discard 121 of 683 samples, and it would wipe out three entire classes, since 2-4-d-injury, cyst-nematode and herbicide-injury would disappear completely and phytophthora-rot would lose 68 of its 88 samples. A model that never sees those diseases cannot predict them.

3. Treat “missing” as its own category rather than imputing. All the predictors are categorical, so a missing value can be coded as an extra level, which is what the code above does. This is the better choice here for two reasons. First, missingness is almost perfectly tied to the class, so the fact that a value is missing is itself useful information, and imputing it would erase that signal. Second, imputation would be unreliable for these samples: a class such as 2-4-d-injury is missing about 80% of its values, so there is very little left on which to base a prediction of the missing entries. As the chapter explains, imputation is a model within a model, and it needs informative predictors to work from. Tree-based models, which Sect. 3.4 notes can handle missing data directly, and naive Bayes both deal with an extra category naturally.

If a method is needed that requires complete numeric data (a neural network or a support vector machine, for example), then K-nearest neighbour or bagged-tree imputation through caret::preProcess, or multiple imputation through the mice package, would be the alternative. In that case I would also add a binary “was missing” indicator for each affected predictor so the information carried by the missingness is kept, and, following the chapter’s advice, I would carry out the imputation inside the resampling loop so that the performance estimates stay honest.