Applied Predictive Modeling: Exercises 3.1 and 3.2
Jonathan Cruz
October 04, 2026
Packages
library(mlbench) # contains the Glass and Soybean datasets
library(ggplot2) # plotting
library(tidyr) # pivot_longer() to reshape data for when we plot
library(dplyr) # data manipulation
library(e1071) # skewness()
library(corrplot) # correlation matrix plot
library(caret) # preProcess() and nearZeroVar() from Chapter 3Exercise 3.1: Glass
The UC Irvine Glass dataset contains 214 glass samples. Each sample has 9 predictors which are the index, the percentage of 8 chemical elements and a class label type that has 6 diffirent types of glass
data(Glass)
str(Glass) # 9 numeric predictors and the factor outcome Type## '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 ...table(Glass$Type) # class counts, types 1, 2 are dominant while 3, 5, 6 are small##
## 1 2 3 5 6 7
## 70 76 17 13 9 29glass_x <- Glass[, 1:9] # keep only the predictors for the analysis below(a) Distributions and relationships between predictors
Distributions. A histogram of each predictor shows
its shape. We want each panel to have their own axis because elements
are measured in diffirent ranges therefore scales = "free"
will accomplish this for us.
# Reshape to long format, one row per predictor and value so ggplot can facet
glass_long <- pivot_longer(glass_x, everything(),
names_to = "predictor", values_to = "value")
ggplot(glass_long, aes(x = value)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white") +
facet_wrap(~ predictor, scales = "free") +
labs(title = "Distribution of each Glass predictor", x = NULL, y = "Count")
Relationships between predictors. A correlation matrix gives us a summary of how every predictors when paired moves together.
glass_cor <- cor(glass_x) # Pearson correlation for every pair
# order = "hclust" groups predictors that are correlated with each other
corrplot(glass_cor, method = "number", type = "upper", order = "hclust",
tl.col = "black", number.cex = 0.8)
The strongest relationship is RI and Ca at 0.81, so we look at that pair directly.
ggplot(Glass, aes(x = Ca, y = RI)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE, color = "red") + # straight-line trend
labs(title = "RI rises with Ca (r = 0.81)")
Answer (a): Na, Al and Si are close to bell-shaped with Al having a slight right tail. RI, Ca and especially K have long right tails. Ba and Fe are zero for most samples, with a few non-zero values. Mg is the opposite, many samples are near 3.5 and have plus a large spike at 0 but nevertheless is left skewed. Most pairs were weakly correlated except for RI and Ca (r = 0.81) which means calcium tends to bend light more.
(b) Outliers and skewness
Skewness. A value near 0 means symmetric, positive means a long right tail, negative means a long left tail. Values beyond about ±1 are considered noticeably skewed.
# Calculate the skewness for each predictor and sort from highest skewed to least skewed
glass_skew <- sort(sapply(glass_x, skewness), decreasing = TRUE)
round(glass_skew, 2)## K Ba Ca Fe RI Al Na Si Mg
## 6.46 3.37 2.02 1.73 1.60 0.89 0.45 -0.72 -1.14Zeros. Counting zeros explains the odd shapes of Ba, Fe, Mg and K.
# How many of the 214 samples in the data are at exactly 0 for each predictor?
sapply(glass_x, function(v) sum(v == 0))## RI Na Mg Al Si K Ca Ba Fe
## 0 0 42 0 0 30 0 176 144Outliers. Boxplots marks any point that’s more than 1.5 × IQR which is Q3 - Q1 beyond the box as an outlier. The predictors are centered and scaled first so all of them fit on one common axis that can be compared.
# scale() centers each predictor to mean 0 and scales it to SD 1
glass_scaled_long <- pivot_longer(as.data.frame(scale(glass_x)), everything(),
names_to = "predictor", values_to = "z")
ggplot(glass_scaled_long, aes(x = predictor, y = z)) +
geom_boxplot(fill = "lightblue", outlier.color = "red") +
labs(title = "Standardized predictors (red points = outliers according to the 1.5 x IQR rule)",
x = NULL, y = "Standard deviations from the mean")
# Count outliers per predictor with the same 1.5 x IQR rule the boxplot uses
count_outliers <- function(v) {
q <- quantile(v, c(0.25, 0.75))
iqr <- q[2] - q[1]
sum(v < q[1] - 1.5 * iqr | v > q[2] + 1.5 * iqr)
}
sapply(glass_x, count_outliers)## RI Na Mg Al Si K Ca Ba Fe
## 17 7 0 18 12 7 26 38 12# The most extreme value is the K reading over 6%, when most are under 1%
Glass[which.max(Glass$K), ]## RI Na Mg Al Si K Ca Ba Fe Type
## 172 1.51316 13.02 0 3.04 70.48 6.21 6.96 0 0 5Answer (b):
There’s definitely outliers, every predictor besides MG has point beyond the boxplot whiskers which K being the most extreme one with 6.21 and is 9 SD above the mean, we have severy skewed predictors K (6.46), Ba (3.37), Ca (2.02), Fe (1.73) and RI (1.60) are right skewed then Mg (-1.14) is left skewed because of that spike at 0. Na, AL and Si are very close to being symmetric
- Caution before deleting outliers: these may be real glass samples, and with only 214 samples (some classes have 9 to 17), removing points could throw away the rare classes we most need.
(c) Transformations that might help
Chapter 3 suggests a Box-Cox transformation for skewed predictors,
however Box-Cox only works on strictly positive values, and Mg, K, Ba
and Fe contain zeros and caret’s preProcess() skips those
columns silently, so in this case Box-Cox would not touch the worst
predictor (K)
# Box-Cox can only be estimated for the positive predictors
bc <- preProcess(glass_x, method = "BoxCox")
bc$method$BoxCox # K, Mg, Ba and Fe are NOT in this list## [1] "RI" "Na" "Al" "Si" "Ca"The Yeo-Johnson transformation is the same idea but it does works with zeros, so we use this method instead.
# Estimate a Yeo-Johnson transformation for each predictor
yj <- preProcess(glass_x, method = "YeoJohnson")
glass_yj <- predict(yj, glass_x) # apply it
# Compare skewness before vs after
skew_compare <- data.frame(
before = round(sapply(glass_x, skewness), 2),
after = round(sapply(glass_yj, skewness), 2)
)
skew_compare## before after
## RI 1.60 1.60
## Na 0.45 -0.01
## Mg -1.14 -0.88
## Al 0.89 0.00
## Si -0.72 -0.72
## K 6.46 -0.07
## Ca 2.02 -0.21
## Ba 3.37 3.37
## Fe 1.73 1.73# Histograms after the transformation took place
pivot_longer(glass_yj, everything(), names_to = "predictor", values_to = "value") |>
ggplot(aes(x = value)) +
geom_histogram(bins = 30, fill = "darkseagreen", color = "white") +
facet_wrap(~ predictor, scales = "free") +
labs(title = "Glass predictors after the Yeo-Johnson transformation", x = NULL, y = "Count")
Chapter 3 has steps that helps us build full preprocessing pipeline to address the problems that we found
# 1. YeoJohnson reduces skew and works with zeros unlike Box-Cox
# 2. put every predictor on the same scale
# 3. pull extreme outliers in closer to the data
glass_pp <- preProcess(glass_x, method = c("YeoJohnson", "center", "scale", "spatialSign"))
glass_pp## Created from 214 samples and 9 variables
##
## Pre-processing:
## - centered (9)
## - ignored (0)
## - scaled (9)
## - spatial sign transformation (9)
## - Yeo-Johnson transformation (5)
##
## Lambda estimates for Yeo-Johnson transformation:
## -0.19, 2.24, 0, -0.98, -1.36glass_ready <- predict(glass_pp, glass_x)Answer (c): Yeo-Johnson instead of Box-Cox absolutely helps with the symmetry as K, Na, Al and Ca came much closer to be symmetric and minimally helped Mg and RI. Centering and scaling are also needed for models that rely on distances or penalties. For example KNN classifies a sample based on nearest neighbors, so predictors with a big spread would dominate the distance calculation. RI, whose SD is about 500 times smaller than Mg’s or Ca’s, would be ignored without scaling.Spatial sign reduces the gravity of extreme outliers like K = 6.21, with the benefit of not deleting any data. For models that struggle with correlated independent variables like RI and Ca at r = 0.81, dropping one of them or using PCA can help.
Exercise 3.2: Soybean
The Soybean dataset has 683 plants, 35 categorical predictors describing the plant and its environment, and an outcome class with 19 diseases.
data(Soybean)
dim(Soybean) # 683 rows, 36 columns ## [1] 683 36str(Soybean[, 1:6]) # predictors are factors, some ordered, coded as 0, 1, 2, ...## 'data.frame': 683 obs. of 6 variables:
## $ Class : Factor w/ 19 levels "2-4-d-injury",..: 11 11 11 11 11 11 11 11 11 11 ...
## $ date : Factor w/ 7 levels "0","1","2","3",..: 7 5 4 4 7 6 6 5 7 5 ...
## $ plant.stand: Ord.factor w/ 2 levels "0"<"1": 1 1 1 1 1 1 1 1 1 1 ...
## $ precip : Ord.factor w/ 3 levels "0"<"1"<"2": 3 3 3 3 3 3 3 3 3 3 ...
## $ temp : Ord.factor w/ 3 levels "0"<"1"<"2": 2 2 2 2 2 2 2 2 2 2 ...
## $ hail : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 1 ...soy_x <- Soybean[, -1] # predictors only(a) Frequency distributions and degenerate predictors
A bar chart of every predictor shows how the samples are spread
across each category. Grey bars (NA) are missing values,
which we examine in part (b).
# Convert every factor to character type so all columns can be stacked into one long table
soy_long <- soy_x |>
mutate(across(everything(), as.character)) |>
pivot_longer(everything(), names_to = "predictor", values_to = "level")
ggplot(soy_long, aes(x = level)) +
geom_bar(fill = "steelblue") +
facet_wrap(~ predictor, scales = "free", ncol = 6) +
labs(title = "Frequency of each category for all 35 Soybean predictors",
x = "Category (NA = missing)", y = "Count") +
theme(strip.text = element_text(size = 7), axis.text = element_text(size = 6))
Chapter 3 defines a degenerate which is near-zero variance predictor as one where both the percentage of unique value is under 10% and the frequency ratio (most common / second most common) is over 20
Every predictor here has only a few categories, so rule 1 is met.
Rule 2 is the deciding factor and nearZeroVar() checks both
rules.
# saveMetrics = TRUE returns the two statistics for each predictor
nzv_stats <- nearZeroVar(soy_x, saveMetrics = TRUE)
# Sort by frequency ratio and show the top 8
head(nzv_stats[order(-nzv_stats$freqRatio), ], 8)## freqRatio percentUnique zeroVar nzv
## mycelium 106.500000 0.2928258 FALSE TRUE
## sclerotia 31.250000 0.2928258 FALSE TRUE
## leaf.mild 26.750000 0.4392387 FALSE TRUE
## shriveling 14.184211 0.2928258 FALSE FALSE
## int.discolor 13.204545 0.4392387 FALSE FALSE
## lodging 12.380952 0.2928258 FALSE FALSE
## leaf.malf 12.311111 0.2928258 FALSE FALSE
## seed.size 9.016949 0.2928258 FALSE FALSE# Just the names of the degenerate predictors
nzv_names <- nearZeroVar(soy_x, names = TRUE)
nzv_names## [1] "leaf.mild" "mycelium" "sclerotia"# Look at the actual counts for the degenerate predictors, using useNA to show missing values too
lapply(soy_x[nzv_names], table, useNA = "ifany")## $leaf.mild
##
## 0 1 2 <NA>
## 535 20 20 108
##
## $mycelium
##
## 0 1 <NA>
## 639 6 38
##
## $sclerotia
##
## 0 1 <NA>
## 625 20 38Answer (a): Three predictors are degenerate which
are mycelium, sclerotia and
leaf.mild, none of them have zero variance as in each
predictor has at least two categories. Many of the other are unbalanced
but not necesarilly degenerate as they are below the 20 cutoff so we
keep those.
(b) Patterns in the missing data
Which predictors are missing most often?
# Percent of rows missing for each predictor
pct_missing <- sort(colMeans(is.na(soy_x)) * 100, decreasing = TRUE)
round(pct_missing, 1)## hail sever seed.tmt lodging germ
## 17.7 17.7 17.7 17.7 16.4
## leaf.mild fruiting.bodies fruit.spots seed.discolor shriveling
## 15.8 15.5 15.5 15.5 15.5
## leaf.shread seed mold.growth seed.size leaf.halo
## 14.6 13.5 13.5 13.5 12.3
## leaf.marg leaf.size leaf.malf fruit.pods precip
## 12.3 12.3 12.3 12.3 5.6
## stem.cankers canker.lesion ext.decay mycelium int.discolor
## 5.6 5.6 5.6 5.6 5.6
## sclerotia plant.stand roots temp crop.hist
## 5.6 5.3 4.5 4.4 2.3
## plant.growth stem date area.dam leaves
## 2.3 2.3 0.1 0.1 0.0data.frame(predictor = names(pct_missing), pct = pct_missing) |>
ggplot(aes(x = reorder(predictor, pct), y = pct)) +
geom_col(fill = "tomato") +
coord_flip() +
labs(title = "Percent missing by predictor", x = NULL, y = "% of plants missing")
~18% comes from a share of plant rows with at lease one missing value, a share of individual cells is much lower.
incomplete <- !complete.cases(Soybean) # TRUE if the row has any NA
sum(incomplete) # 121 plants have at least one missing value## [1] 121round(mean(incomplete) * 100, 1) # 17.7% of rows which is the book's "roughly 18%"## [1] 17.7round(mean(is.na(soy_x)) * 100, 1) # 9.5% of all individual cells## [1] 9.8Check for is missingness related to the class? Let’s Count complete vs incomplete plants within each disease.
# For each class number of plants, number with any missing, and percent
missing_by_class <- Soybean |>
mutate(incomplete = !complete.cases(Soybean)) |>
group_by(Class) |>
summarise(plants = n(),
with_missing = sum(incomplete),
pct_with_missing = round(mean(incomplete) * 100, 1)) |>
arrange(desc(pct_with_missing))
missing_by_class## # A tibble: 19 × 4
## Class plants with_missing pct_with_missing
## <fct> <int> <int> <dbl>
## 1 2-4-d-injury 16 16 100
## 2 cyst-nematode 14 14 100
## 3 diaporthe-pod-&-stem-blight 15 15 100
## 4 herbicide-injury 8 8 100
## 5 phytophthora-rot 88 68 77.3
## 6 alternarialeaf-spot 91 0 0
## 7 anthracnose 44 0 0
## 8 bacterial-blight 20 0 0
## 9 bacterial-pustule 20 0 0
## 10 brown-spot 92 0 0
## 11 brown-stem-rot 44 0 0
## 12 charcoal-rot 20 0 0
## 13 diaporthe-stem-canker 20 0 0
## 14 downy-mildew 20 0 0
## 15 frog-eye-leaf-spot 91 0 0
## 16 phyllosticta-leaf-spot 20 0 0
## 17 powdery-mildew 20 0 0
## 18 purple-seed-stain 20 0 0
## 19 rhizoctonia-root-rot 20 0 0A heatmap shows the full pattern, the percent of each predictor that is missing within each class.
# For each class, the percent missing in each predictor
heat_data <- Soybean |>
group_by(Class) |>
summarise(across(everything(), ~ mean(is.na(.x)) * 100)) |>
pivot_longer(-Class, names_to = "predictor", values_to = "pct_missing")
ggplot(heat_data, aes(x = predictor, y = Class, fill = pct_missing)) +
geom_tile(color = "white") +
scale_fill_gradient(low = "white", high = "darkred", name = "% missing") +
labs(title = "Missing data is concentrated in 5 of the 19 classes", x = NULL, y = NULL) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 9),
axis.text.y = element_text(size = 9))
Answer (b): There are some predictors that are more
likely to be missing, hail, sever,
seed.tmt and lodging are each missing for
17.7% of plants, followed by germ (16.4%) and
leaf.mild (15.8%). Predictors go missing in groups like the
seed and fruit measurements (seed, seed.size,
fruit.spots, shriveling…) tend to be missing
together, as do the leaf measurements. Patterns are strongly tied to the
class as all incomplete pants belongs to just 5 of the 19 classes and in
4 of them which are 2-4-d-injury,
cyst-nematode, diaporthe-pod-&-stem-blight
and herbicide-injury every plant has missing values.The
other 14 classes have no missing values at all.
(c) Strategy for handling missing data
First, compare the simple options so the choice is justified by numbers not guessing.
# 1, drop every row with a missing value
rows_complete <- Soybean[complete.cases(Soybean), ]
# 2. drop every predictor with any missing value
# drop = FALSE keeps it a data frame even if only one column survives
cols_complete <- soy_x[, colSums(is.na(soy_x)) == 0, drop = FALSE]
names(cols_complete) # the only predictors with no missing values## [1] "leaves"data.frame(
option = c("Keep everything (no fix)",
"1. Drop rows with any NA",
"2. Drop predictors with any NA"),
plants = c(nrow(Soybean), nrow(rows_complete), nrow(cols_complete)),
predictors = c(ncol(soy_x), ncol(rows_complete) - 1, ncol(cols_complete)),
classes = c(nlevels(Soybean$Class),
n_distinct(rows_complete$Class), # how many diseases survive?
nlevels(Soybean$Class))
)## option plants predictors classes
## 1 Keep everything (no fix) 683 35 19
## 2 1. Drop rows with any NA 562 35 15
## 3 2. Drop predictors with any NA 683 1 19Neither deletion options work as dropping rows loses 121 plants and
completely wipes out 4 diseases and the model would not be able to drop
them. Dropping predictors leaves us with only 1 out of the 35
predictors(leaves) which basically throws away almost all
of the information
Chosen strategy:
Remove the 3 degenerate predictors discussed in part a because they basically hold no information and then treat “missing” as its own category. The predictors are already categorical and this would let us keep all 683 plants and all 19 classes. It also allows for the model to use missingness as information.
# 1st step -> remove the near zero variance predictors found in part a
soy_clean <- Soybean[, !(names(Soybean) %in% nzv_names)]
# 2nd step -> helper that turns NA into an explicit "missing" category
add_missing_level <- function(f) {
f <- factor(f, ordered = FALSE) # drop the ordering so we can add a level
levels(f) <- c(levels(f), "missing") # add the new category
f[is.na(f)] <- "missing" # fill NA with it
f
}
# Apply the helper to every predictor
soy_clean[, -1] <- lapply(soy_clean[, -1], add_missing_level)
dim(soy_clean) # 683 plants, 33 columns ## [1] 683 33sum(is.na(soy_clean)) # 0 missing values remain## [1] 0nlevels(soy_clean$Class) # all 19 classes kept## [1] 19table(soy_clean$hail) # "missing" is now a real category##
## 0 1 missing
## 435 127 121Answer (c): My strategy has two steps which first I
removed the three near zero variance predictors (leaf.mild,
mycelium and sclerotia), since they barely
vary and essentially adds almost no information whatsoever. Second,
instead of deleting rows or filling in values, I categorized every
missing value as its own “missing” category. This kept all 683 plants
and all 19 classes and because all the missing values come entirely from
five diseases, the model can use that as a recognizable pattern . If the
missing values had been random instead, imputing them with KNN or bagged
trees would have been a reasonable alternative.