Both exercises are from Chapter 3 of Kuhn and Johnson, Applied Predictive Modeling, and use data sets from the mlbench package: Glass for 3.1 and Soybean for 3.2.

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, including the refractive index and 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 ...

The structure output shows six levels in Type, not seven. Type 4 has no samples in this copy of the data, so the classification problem here has six classes.

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

# one histogram per predictor in a 3 x 3 grid
# Type (column 10) is the outcome, so it is left out
par(mfrow = c(3, 3))
for (v in names(Glass)[1:9]) {
  hist(Glass[[v]], breaks = 25, main = v, xlab = "", col = "lightblue")
}

par(mfrow = c(1, 1))

Reading the histograms one at a time:

  • RI, Al, Ca, Na, Si each have a single main hump, with a tail on one side for most of them.
  • Mg is different. There is a tall group of samples at exactly zero and a second group between about 3 and 4.5, with little in between, so it is bimodal.
  • K, Ba and Fe are piled up at or near zero with a few much larger values stretching to the right.

The class counts matter for the later modeling too:

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

Types 1 and 2 hold most of the samples, and Type 6 is the smallest class.

Now the relationships between predictors:

# correlations among the nine predictors; hclust order groups similar ones
glass_cor <- cor(Glass[, 1:9])
corrplot(glass_cor, order = "hclust")

round(glass_cor, 2)
##       RI    Na    Mg    Al    Si     K    Ca    Ba    Fe
## RI  1.00 -0.19 -0.12 -0.41 -0.54 -0.29  0.81  0.00  0.14
## Na -0.19  1.00 -0.27  0.16 -0.07 -0.27 -0.28  0.33 -0.24
## Mg -0.12 -0.27  1.00 -0.48 -0.17  0.01 -0.44 -0.49  0.08
## Al -0.41  0.16 -0.48  1.00 -0.01  0.33 -0.26  0.48 -0.07
## Si -0.54 -0.07 -0.17 -0.01  1.00 -0.19 -0.21 -0.10 -0.09
## K  -0.29 -0.27  0.01  0.33 -0.19  1.00 -0.32 -0.04 -0.01
## Ca  0.81 -0.28 -0.44 -0.26 -0.21 -0.32  1.00 -0.11  0.12
## Ba  0.00  0.33 -0.49  0.48 -0.10 -0.04 -0.11  1.00 -0.06
## Fe  0.14 -0.24  0.08 -0.07 -0.09 -0.01  0.12 -0.06  1.00

Most pairs are weak (the faint circles). A handful are moderate: RI with Si (negative), Al with Ba (positive), and Mg with Al, Ca and Ba (all negative). The one strong relationship is RI and Ca, which are strongly positively related. A scatterplot confirms it:

plot(Glass$Ca, Glass$RI, xlab = "Ca", ylab = "RI", pch = 16,
     col = "steelblue")

RI and Ca rise together, so for a model that dislikes correlated predictors (the chapter’s example is linear regression) one of the two carries mostly redundant information. The other pairwise correlations are smaller, though that alone does not rule out relationships involving several predictors at once.

To see whether any predictor separates the classes, here are boxplots by Type:

par(mfrow = c(3, 3))
for (v in names(Glass)[1:9]) {
  boxplot(Glass[[v]] ~ Glass$Type, main = v, xlab = "Type", ylab = "",
          col = "lightblue")
}

par(mfrow = c(1, 1))

Mg is low for types 5, 6 and 7 and high for the others. Ba has a clearly positive median only for type 7, while the other types have a median of zero (a few of those types still contain positive Ba values, which show as dots). Al, Na and Ca have centers that shift for some types (Ca for type 5, Na for types 6 and 7), so those look informative for classification. Si and Fe overlap heavily across all six types.

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

par(mfrow = c(3, 3))
for (v in names(Glass)[1:9]) {
  boxplot(Glass[[v]], main = v, col = "lightblue")
}

par(mfrow = c(1, 1))
# skewness(): near 0 is symmetric, large positive is right-skewed
# boxplot.stats()$out counts points beyond 1.5 x IQR (the boxplot rule)
glass_summary <- data.frame(
  skewness = sapply(Glass[, 1:9], skewness),
  boxplot_outliers = sapply(Glass[, 1:9],
                            function(x) length(boxplot.stats(x)$out)),
  zeros = sapply(Glass[, 1:9], function(x) sum(x == 0))
)
round(glass_summary, 2)
##    skewness boxplot_outliers zeros
## RI     1.60               17     0
## Na     0.45                7     0
## Mg    -1.14                0    42
## Al     0.89               18     0
## Si    -0.72               12     0
## K      6.46                7    30
## Ca     2.02               26     0
## Ba     3.37               38   176
## Fe     1.73               12   144

Outliers. Yes. The boxplots mark points beyond the whiskers for every predictor except Mg, and the table gives the counts. K, Ba, Ca and RI stand out most: K has two samples at its largest value, far above the rest, and Ba has a long run of points above a box that is squeezed flat at zero. Because 176 of the 214 Ba values are zero, the box has no height, so every positive Ba value gets flagged. That is how the boxplot rule works here, not evidence that each positive value is a mistake. Following the chapter’s caution, I would not delete any of these on sight. With only 214 samples, a high Ba or K value may belong to a real, rare kind of glass, and most of the positive Ba values do belong to Type 7 (table below). They should be checked as possible recording errors, and otherwise kept.

# the samples at the largest K value, and the count of positive Ba values
# in each type
Glass[Glass$K == max(Glass$K), c("K", "Type")]
##        K Type
## 172 6.21    5
## 173 6.21    5
with(Glass, table(Type, Ba_positive = Ba > 0))
##     Ba_positive
## Type FALSE TRUE
##    1    67    3
##    2    70    6
##    3    16    1
##    5    11    2
##    6     9    0
##    7     3   26

Skewness. Yes. K has by far the largest skewness, followed by Ba, Ca, Fe and RI, all strongly right-skewed, with Al moderately so. Mg and Si are skewed to the left. Na is only mildly skewed. Part of the skew for Ba, Fe and K comes from the zeros in the zeros column: both Ba and Fe are zero for most samples, with a tail of positive values, so a large share of their distribution is a single repeated value.

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

The chapter’s remedy for right skew is a log, square root or inverse transformation, with Box-Cox choosing the power from the data. Box-Cox requires strictly positive values, so it can only be applied to the predictors with no zeros.

# BoxCoxTrans() estimates lambda by maximum likelihood
# it cannot handle zeros, so only the predictors with none are used
positive_cols <- c("RI", "Na", "Al", "Si", "Ca")
sapply(Glass[, positive_cols], function(x) BoxCoxTrans(x)$lambda)
##   RI   Na   Al   Si   Ca 
## -2.0 -0.1  0.5  2.0 -1.1

Al and Ca are two promising candidates: both are right-skewed, and their estimated lambdas are close to a square root (0.5) and an inverse (-1), respectively. Si’s estimate of 2 and RI’s estimate of -2 sit at the ends of the range that was searched, so I focus on Al and Ca here. Here is the effect for those two:

al_bc <- BoxCoxTrans(Glass$Al)
ca_bc <- BoxCoxTrans(Glass$Ca)
al_new <- predict(al_bc, Glass$Al)
ca_new <- predict(ca_bc, Glass$Ca)

par(mfrow = c(2, 2))
hist(Glass$Al, breaks = 25, main = "Al, original", xlab = "",
     col = "lightblue")
hist(al_new, breaks = 25, main = "Al, Box-Cox", xlab = "",
     col = "lightgreen")
hist(Glass$Ca, breaks = 25, main = "Ca, original", xlab = "",
     col = "lightblue")
hist(ca_new, breaks = 25, main = "Ca, Box-Cox", xlab = "",
     col = "lightgreen")

par(mfrow = c(1, 1))

c(Al_before = skewness(Glass$Al), Al_after = skewness(al_new),
  Ca_before = skewness(Glass$Ca), Ca_after = skewness(ca_new))
##   Al_before    Al_after   Ca_before    Ca_after 
##  0.89461042  0.09105899  2.01844629 -0.19395573

The transformed histograms are noticeably more symmetric and the skewness statistics fall toward zero, so these transformations reduce the skewness. Some extreme values remain in the tails. This shows the skew going down; no classification model was fitted, so it does not show better predictions.

For the predictors Box-Cox cannot touch, the choices are different:

  • K, Ba, Fe. These contain zeros, so a log and the Box-Cox estimation fail. A square root is defined at zero and is one of the chapter’s simple candidates; the table below shows what it does to their skewness. For Ba and Fe the larger issue is that most values are zero, which no power transformation removes. Recoding them as present/absent would throw away the measured amounts (the chapter cautions against binning for that reason), so I would keep the numeric values and compare an indicator version only if a model needed it. A tree splits on thresholds and is not bothered by skew or zeros.
  • Mg. It has a group at zero, another around 3 to 4.5, and negative skewness. A power transformation does not remove the repeated zeros, so making Mg symmetric is not a sufficient goal on its own, and a tree or other flexible model can use the zero group directly.
  • Outliers. If the chosen model is sensitive to them (a linear boundary, for example), the spatial sign transformation from the chapter pulls extreme samples toward the rest. The chapter says to center and scale first, since it works on all predictors together.
  • Centering and scaling. The predictors have very different spreads and units (the histograms show RI covering a tiny range and Si a much larger one), so centering and scaling puts them on a common standard-deviation scale for distance-based methods and PCA.
  • Correlation. RI and Ca are the one strongly correlated pair (about 0.81). For a correlation-sensitive model I would consider dropping one of them or replacing the two with a principal component, after centering and scaling.
# skewness before and after a square root, for the predictors with zeros
zero_cols <- Glass[, c("K", "Ba", "Fe")]
round(rbind(original = sapply(zero_cols, skewness),
            sqrt = sapply(zero_cols, function(x) skewness(sqrt(x)))), 2)
##             K   Ba   Fe
## original 6.46 3.37 1.73
## sqrt     0.86 2.34 1.04

The square root lowers the skewness for all three, though Ba stays clearly skewed because of its many zeros.

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 (e.g., temperature, precipitation) and plant conditions (e.g., left spots, mold growth). The outcome labels consist of 19 distinct classes.

data(Soybean)
# ?Soybean for details
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?

# bar chart of level counts for each of the 35 predictors
# (column 1 is the outcome, Class)
par(mfrow = c(7, 5), mar = c(2.5, 2.5, 2, 0.5))
for (v in names(Soybean)[-1]) {
  counts <- table(Soybean[[v]], useNA = "ifany")
  labels <- names(counts)
  labels[is.na(labels)] <- "NA"
  barplot(counts, names.arg = labels, main = v, col = "lightblue",
          cex.names = 0.8)
}

par(mfrow = c(1, 1))

In many panels one level is much taller than the rest, and some panels have an extra bar labeled NA for the missing values. The chapter’s definition of a degenerate predictor is a zero-variance or near-zero-variance one: very few unique values, and one value overwhelmingly more common than the next. nearZeroVar() applies the chapter’s two rules (a low fraction of unique values, about 10% or less, and a ratio of the most common to the second most common value of about 20 or more; caret’s default cutoff is 95/5 = 19). The outcome Class is left out of the check:

nzv_info <- nearZeroVar(Soybean[, -1], saveMetrics = TRUE)
nzv_info[nzv_info$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

Three predictors are flagged: leaf.mild, mycelium and sclerotia. None has zero variance (every one has at least two levels), but each is dominated by a single level. mycelium is the most extreme: nearly every sample has level 0, and only a handful have level 1.

table(Soybean$leaf.mild, useNA = "ifany")
## 
##    0    1    2 <NA> 
##  535   20   20  108
table(Soybean$mycelium, useNA = "ifany")
## 
##    0    1 <NA> 
##  639    6   38
table(Soybean$sclerotia, useNA = "ifany")
## 
##    0    1 <NA> 
##  625   20   38

These predictors have highly imbalanced levels. The chapter warns that resampling can make a rare level disappear from a resample, which causes problems for some models. A tree cannot split on a constant predictor, but it may still use a near-zero-variance one if its rare values help separate the classes, so whether to drop them depends on the model.

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

Based on (a) and (b), my strategy has three parts.

  1. Keep every sample. Removing incomplete rows would wipe out whole classes, as shown above.
  2. Do not drop predictors for missingness. The worst predictors are missing for under a fifth of samples, and the missingness itself separates the classes, so removing those columns would throw away signal.
  3. Treat “missing” as its own level. Section 3.4 emphasizes that missingness can be informative, and Section 3.6 (Table 3.2) shows an “Unknown” category kept for a categorical predictor. I propose replacing NA with an explicit “missing” category for these predictors. This preserves which measurements are absent while keeping every sample and its observed values. It does not estimate the unobserved value, and nothing here shows it predicts better than another treatment, since that would need a model and resampling.
# addNA() adds NA as a real factor level
# the level is then renamed so it prints as "missing"
soy_clean <- Soybean
for (v in names(soy_clean)[-1]) {
  soy_clean[[v]] <- addNA(soy_clean[[v]])
  levels(soy_clean[[v]])[is.na(levels(soy_clean[[v]]))] <- "missing"
}
# check: no NA values remain in the predictors
sum(is.na(soy_clean[, -1]))
## [1] 0

The near-zero variance check has to be re-run on the cleaned data, because “missing” is now a level and changes the frequency ratios:

nzv_clean <- nearZeroVar(soy_clean[, -1], names = TRUE)
nzv_clean
## character(0)
soy_final <- soy_clean[, !(names(soy_clean) %in% nzv_clean)]
dim(soy_final)
## [1] 683  36

No predictor meets the near-zero-variance criteria once missingness is its own category, so soy_final keeps all 35 predictors plus Class. The explicit category does not create new measurements, and it does not show that every predictor is useful.

The three predictors flagged earlier drop below the cutoff because the missing values now count as a level:

# the three predictors flagged earlier, checked on the original data
# and again after "missing" became its own level
flagged <- c("leaf.mild", "mycelium", "sclerotia")
keep <- c("freqRatio", "percentUnique")
before <- nearZeroVar(Soybean[, -1], saveMetrics = TRUE)[flagged, keep]
after <- nearZeroVar(soy_clean[, -1], saveMetrics = TRUE)[flagged, keep]
out <- data.frame(before, after)
names(out) <- c("ratio_before", "unique_pct_before",
                "ratio_after", "unique_pct_after")
round(out, 2)
##           ratio_before unique_pct_before ratio_after unique_pct_after
## leaf.mild        26.75              0.44        4.95             0.59
## mycelium        106.50              0.29       16.82             0.44
## sclerotia        31.25              0.29       16.45             0.44

All three ratios fall below the cutoff of about 19 once missing values are counted as a level, because the missing group is now the second most common value. These predictors are still dominated by one level, and mycelium and sclerotia sit close to the cutoff.

Section 3.4 also describes estimating a missing predictor from the other predictors, for example with K-nearest neighbors. For categorical variables like these, any such imputation has to return valid category levels, since averaging numeric category codes is not appropriate. It should be fit inside each training resample, using the predictors and not the disease label, so that it could be applied to new samples and the performance estimate stays honest. Models that need numeric inputs could use dummy variables for the explicit categories, including “missing” (Section 3.6). Some tree implementations handle missing predictors directly, which is another option.