Rpubs URL: https://rpubs.com/umaisabdullah/Data624_Homework4

Objective

Complete and submit Exercises 3.1 and 3.2 from Applied Predictive Modeling by Kuhn and Johnson. We have to submit both our Rpubs link as well as attach the .pdf file with our code.

Overview

The previous homeworks were about time series, where the order of the observations is what carries the information. These two exercises are a different kind of problem. The samples have no time order at all, so instead of looking for trend and seasonality I am looking at the shape of each predictor and at the problems that would cause trouble for a model later on.

The two exercises cover different problems. Exercise 3.1 uses the Glass data, where the predictors are numbers and the issues are skewness, outliers and correlation. Exercise 3.2 uses the Soybean data, where the predictors are categories and the issues are predictors that barely vary and a large amount of missing data.

What the two have in common is the order of work. In both cases I look at the data first, put numbers on whatever looks wrong, and only then decide what to do about it. Choosing the fix before understanding the problem is how you end up applying a transformation that does nothing.

Pre-Requisite

Both datasets come from the mlbench package, so nothing needs to be downloaded into the working directory. The packages used are:

mlbench contains the Glass and Soybean datasets.
e1071 provides the skewness() function used to measure how lopsided each predictor is.
corrplot draws the correlation matrix.
caret provides nearZeroVar(), which is used to find degenerate predictors in Exercise 3.2.
ggplot2, tidyr and dplyr are used for the plots and data reshaping.

Exercise 3.1

Question. 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. The data can be accessed via:

library(mlbench)
data(Glass)
str(Glass)

1. Using visualizations, explore the predictor variables to understand their distributions as well as the relationships between predictors.
1. Do there appear to be any outliers in the data? Are any predictors skewed?
1. Are there any relevant transformations of one or more predictors that might improve the classification model?

My thought process. The three parts build on each other, so I work through them in order rather than jumping straight to the transformations. Part a tells me the shape of each predictor. Part b puts numbers on the problems I can see in those shapes. Part c is where I decide what to actually do about them.

The one thing I want to avoid is transforming everything out of habit. A transformation should fix a specific problem, so I need to know what the problems are first. That is why part c comes last.

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 ...

# How many samples are in each class?
table(Glass$Type)

## 
##  1  2  3  5  6  7 
## 70 76 17 13  9 29

Before looking at the predictors, the class counts are worth noting. The classes are very unbalanced. Types 1 and 2 have 70 and 76 samples, while Type 6 has only 9. That is not what the question asks about, but it is important context, because a classifier can score well on the big classes and still be useless on the small ones.

Note also that the question says seven categories but the factor has six levels. There is no Type 4 in the data, so only six classes actually appear.

(a) Distributions and relationships between predictors

My thought process. There are nine predictors, so drawing nine separate plots one at a time would be slow to read. Instead I reshape the data into long form and use facets, which puts every predictor on one page and makes comparison easy. I use free scales because the predictors are on very different ranges. Silicon sits around 72 percent while iron is near zero, so forcing them onto a shared axis would flatten most of the panels into nothing.

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

ggplot(glass_long, aes(x = Value)) +
  geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
  facet_wrap(~ Predictor, scales = "free", ncol = 3) +
  labs(title = "Distribution of each predictor",
       x = "Value", y = "Count")

What the histograms show. The nine predictors fall into three clear groups, and the grouping matters more than any single panel.

The first group is roughly symmetric. Na, Al and Si look reasonably bell shaped, though Si leans slightly to the left. These are the predictors that need the least attention.

The second group is strongly right skewed. K, Ba, Fe and Ca all pile up at the low end with a long tail stretching to the right. K is the most extreme, with almost everything near zero and a couple of values far out on the right.

The third group is the unusual one. Mg does not have a single peak. It has a large clump of values near 3.5 and a second clump sitting exactly at zero, with a gap in between. Ba and Fe do something similar, where most samples are exactly zero and the rest spread out above it.

That third group is worth pausing on, because a spike at exactly zero usually does not mean a small measurement. It means the element was not detected at all, which is a different kind of thing than a low reading.

# How many samples sit at exactly zero for each predictor?
Glass |>
  select(-Type) |>
  summarise(across(everything(), \(x) sum(x == 0))) |>
  pivot_longer(everything(), names_to = "Predictor", values_to = "Zeros") |>
  arrange(desc(Zeros))

## # A tibble: 9 × 2
##   Predictor Zeros
##   <chr>     <int>
## 1 Ba          176
## 2 Fe          144
## 3 Mg           42
## 4 K            30
## 5 RI            0
## 6 Na            0
## 7 Al            0
## 8 Si            0
## 9 Ca            0

This confirms what the histograms hinted at. Ba is zero in 176 of the 214 samples, Fe in 144, Mg in 42, and K in 30. So four of the nine predictors are mostly zero. These are not ordinary continuous measurements, and that fact drives most of my answer to part c.

ggplot(glass_long, aes(y = Value)) +
  geom_boxplot(fill = "lightsteelblue") +
  facet_wrap(~ Predictor, scales = "free_y", ncol = 3) +
  labs(title = "Boxplot of each predictor", y = "Value") +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

The boxplots say the same thing as the histograms but make the extreme values easier to count, which is what part b asks about.

Relationships between predictors

cor_matrix <- cor(Glass[, 1:9])

corrplot(cor_matrix, method = "color", type = "upper",
         addCoef.col = "black", number.cex = 0.7,
         tl.col = "black", tl.srt = 45,
         title = "Correlation between predictors",
         mar = c(0, 0, 2, 0))

# List the strongest pairs so I am reading numbers, not colours.
cm <- cor(Glass[, 1:9])
cm[lower.tri(cm, diag = TRUE)] <- NA

as.data.frame(as.table(cm)) |>
  filter(!is.na(Freq), abs(Freq) > 0.4) |>
  arrange(desc(abs(Freq))) |>
  rename(Predictor_1 = Var1, Predictor_2 = Var2, Correlation = Freq) |>
  mutate(Correlation = round(Correlation, 3))

##   Predictor_1 Predictor_2 Correlation
## 1          RI          Ca       0.810
## 2          RI          Si      -0.542
## 3          Mg          Ba      -0.492
## 4          Mg          Al      -0.482
## 5          Al          Ba       0.479
## 6          Mg          Ca      -0.444
## 7          RI          Al      -0.407

What the correlations show. Most pairs are weakly related, but a few stand out.

By far the strongest is RI and Ca at 0.81. That one relationship explains about 66 percent of the variation in either variable, so these two predictors carry a lot of the same information. This makes physical sense. The refractive index of glass depends on what it is made of, and calcium is one of the main ingredients, so the two move together.

After that, RI and Si are negatively related at -0.54, and there is a cluster of moderate relationships involving Mg, which is negatively related to Ba, Al and Ca.

ggplot(Glass, aes(x = Ca, y = RI)) +
  geom_point(alpha = 0.6, colour = "steelblue") +
  geom_smooth(method = "lm", se = FALSE, colour = "darkred") +
  labs(title = "Refractive index against calcium",
       subtitle = "The strongest relationship in the data, correlation 0.81",
       x = "Ca (percent)", y = "Refractive index")

The scatterplot confirms it is a genuine straight line relationship and not an artifact of a few stray points.

ggplot(Glass, aes(x = Type, y = Mg, fill = Type)) +
  geom_boxplot(show.legend = FALSE) +
  labs(title = "Magnesium by glass type",
       x = "Glass type", y = "Mg (percent)")

I also looked at the predictors against the class, since the whole point is classification. Magnesium separates the types quite well. Types 1 and 3 sit high, around 3.5, while types 5, 6 and 7 sit much lower. So the odd two clump shape of Mg is not noise. It reflects a real difference between kinds of glass.

(b) Outliers and skewness

My thought process. Both of these questions can be answered by eye from the plots above, but eyes are easy to fool. So I put numbers on both. For skewness I use the skewness() function, and for outliers I use the standard rule of counting points more than 1.5 times the interquartile range beyond the quartiles. Using a stated rule means I am being consistent across the nine predictors instead of judging each one differently.

Skewness

skew_values <- Glass |>
  select(-Type) |>
  summarise(across(everything(), skewness)) |>
  pivot_longer(everything(), names_to = "Predictor", values_to = "Skewness") |>
  arrange(desc(Skewness)) |>
  mutate(Skewness = round(Skewness, 3))

skew_values

## # A tibble: 9 × 2
##   Predictor Skewness
##   <chr>        <dbl>
## 1 K            6.46 
## 2 Ba           3.37 
## 3 Ca           2.02 
## 4 Fe           1.73 
## 5 RI           1.60 
## 6 Al           0.895
## 7 Na           0.448
## 8 Si          -0.72 
## 9 Mg          -1.14

Are any predictors skewed? Yes, several. A skewness near 0 means symmetric. As a rough guide, beyond about 1 in either direction is considered strongly skewed.

By that standard, K is extremely right skewed at 6.46, which is far beyond anything else in the data. Ba at 3.37, Ca at 2.02, Fe at 1.73 and RI at 1.60 are all strongly right skewed as well. Mg is the only meaningfully left skewed predictor at -1.14, which comes from that clump of zeros pulling the tail to the left.

Na at 0.45 and Si at -0.72 are close enough to symmetric that I would leave them alone. Al at 0.895 sits right at the edge of the threshold; it is borderline rather than clearly skewed, so I would not transform it on its own but would revisit it if a later model’s diagnostics suggest it is still causing trouble.

Outliers

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

Glass |>
  select(-Type) |>
  summarise(across(everything(), outlier_count)) |>
  pivot_longer(everything(), names_to = "Predictor", values_to = "Outliers") |>
  arrange(desc(Outliers))

## # A tibble: 9 × 2
##   Predictor Outliers
##   <chr>        <int>
## 1 Ba              38
## 2 Ca              26
## 3 Al              18
## 4 RI              17
## 5 Si              12
## 6 Fe              12
## 7 Na               7
## 8 K                7
## 9 Mg               0

Do there appear to be any outliers? Yes, but they need careful reading. The rule flags 137 values in total across the nine predictors. Ba has the most at 38, followed by Ca at 26, Al at 18 and RI at 17. Mg is the only predictor with none.

Here is the important part. The Ba and Fe counts are misleading. Because those predictors are zero in most samples, the interquartile range is nearly zero, so almost any non zero reading gets flagged as an outlier automatically. Those are not strange measurements. They are ordinary readings from a predictor whose distribution is mostly zeros. Applying a rule without thinking about the shape of the data would give the wrong conclusion here.

The K values are a different story and are worth looking at directly.

# The two largest potassium values.
Glass |>
  mutate(row = row_number()) |>
  filter(K > 3) |>
  select(row, RI, Na, Mg, Al, Si, K, Ca, Ba, Fe, Type)

##     row      RI    Na Mg   Al    Si    K   Ca Ba Fe Type
## 172 172 1.51316 13.02  0 3.04 70.48 6.21 6.96  0  0    5
## 173 173 1.51321 13.00  0 3.02 70.70 6.21 6.93  0  0    5

# How much do those two rows drive the skewness?
c(with_all_data     = round(skewness(Glass$K), 2),
  without_those_two = round(skewness(Glass$K[Glass$K < 3]), 2))

##     with_all_data without_those_two 
##              6.46              1.72

What I found in K. Two samples, rows 172 and 173, both have K equal to 6.21 when the rest of the data sits below 1.6. Those two rows are also nearly identical to each other across every other predictor, which makes me think they are duplicate or near duplicate measurements of the same piece of glass rather than two independent samples.

Those two values alone drive the skewness from 1.72 up to 6.46. So most of the extreme skewness in K comes from two rows, not from the general shape of the predictor. That is worth knowing, because it changes the right fix. A transformation treats the whole column, but the real issue here is two specific rows.

(c) Relevant transformations

My thought process. Now that I know what the problems are, I can match a fix to each one instead of applying the same treatment everywhere. There are three separate problems in this data and they need three different answers, which is the main point I want to make here.

Problem 1: the mostly zero predictors

Ba, Fe, Mg and K are zero in a large share of the samples. A Box-Cox transformation cannot be applied directly to any of them, because Box-Cox needs strictly positive values and taking a log of zero is undefined.

compare_skew <- function(x, name) {
  data.frame(
    Predictor = name,
    Raw       = round(skewness(x), 2),
    Log1p     = round(skewness(log1p(x)), 2),
    SquareRoot = round(skewness(sqrt(x)), 2)
  )
}

bind_rows(
  compare_skew(Glass$K,  "K"),
  compare_skew(Glass$Ba, "Ba"),
  compare_skew(Glass$Fe, "Fe"),
  compare_skew(Glass$Ca, "Ca"),
  compare_skew(Glass$RI, "RI")
)

##   Predictor  Raw Log1p SquareRoot
## 1         K 6.46  1.95       0.86
## 2        Ba 3.37  2.69       2.34
## 3        Fe 1.73  1.56       1.04
## 4        Ca 2.02  1.15       1.55
## 5        RI 1.60  1.60       1.60

The table shows what a shifted log or a square root actually achieves. For K, the square root brings skewness from 6.46 down to 0.86, which is a genuine improvement. For Ba the same transformation only moves it from 3.37 to 2.34, which is not much help.

The reason Ba barely improves is that no transformation can fix its real problem. The issue is not the scale of the values, it is that 176 of 214 samples are exactly zero. Stretching or squashing the axis does not change the fact that most of the mass sits on a single point.

What I would do instead for Ba and Fe. Rather than transforming them, I would turn each one into a yes or no indicator of whether the element was detected at all. That matches what the data is really telling us and it is far more useful for classification.

Glass |>
  group_by(Type) |>
  summarise(n = n(),
            proportion_Ba_present = round(mean(Ba > 0), 3),
            proportion_Mg_zero    = round(mean(Mg == 0), 3))

## # A tibble: 6 × 4
##   Type      n proportion_Ba_present proportion_Mg_zero
##   <fct> <int>                 <dbl>              <dbl>
## 1 1        70                 0.043              0    
## 2 2        76                 0.079              0.118
## 3 3        17                 0.059              0    
## 4 5        13                 0.154              0.538
## 5 6         9                 0                  0.333
## 6 7        29                 0.897              0.793

This is the strongest result in my analysis. Barium is present in 89.7 percent of Type 7 samples and in almost none of the others. So the simple question of whether barium was detected at all is close to a direct marker for Type 7 glass. Magnesium works in a similar way in the opposite direction, with Mg being zero in 79 percent of Type 7 and 54 percent of Type 5, but never in Types 1 or 3.

A model given a plain Ba percentage has to work that out on its own from a column that is mostly zeros. A model given a Ba present indicator gets the useful signal directly. That is a better transformation than any power transformation would be.

Problem 2: the genuinely skewed continuous predictors

Ca and RI are skewed but have no zeros, so a Box-Cox transformation applies to them normally.

# Box-Cox only works on strictly positive predictors.
positive_preds <- c("RI", "Na", "Al", "Si", "Ca")

sapply(positive_preds, function(nm) {
  x  <- Glass[[nm]]
  bc <- MASS::boxcox(x ~ 1, lambda = seq(-3, 3, 0.1), plotit = FALSE)
  round(bc$x[which.max(bc$y)], 2)
})

##   RI   Na   Al   Si   Ca 
## -3.0 -0.1  0.5  3.0 -1.1

data.frame(
  Original = Glass$Ca,
  Logged   = log1p(Glass$Ca)
) |>
  pivot_longer(everything(), names_to = "Version", values_to = "Value") |>
  ggplot(aes(x = Value)) +
  geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
  facet_wrap(~ Version, scales = "free") +
  labs(title = "Calcium before and after a log transformation",
       x = "Value", y = "Count")

For calcium a log transformation cuts the skewness from 2.02 to about 1.15, the same log1p transform used in the comparison table above, and makes the distribution noticeably more even. That is a reasonable improvement and I would keep it.

RI is the odd one out. The suggested lambda runs off to the edge of the range I searched, which is a sign that the method is struggling rather than finding a good answer. The reason is that RI barely varies at all, ranging only from 1.511 to 1.534 with a standard deviation of 0.003. Its apparent skewness comes from a handful of points in a very narrow band. I would centre and scale it rather than transform its shape.

Problem 3: the correlation between RI and Ca

RI and Ca correlate at 0.81, so they carry much of the same information. For models that are sensitive to correlated predictors, such as linear discriminant analysis or logistic regression, I would either drop one of them or replace the pair with principal components.

Tree based models do not need this, since they are not troubled by correlated predictors. So this decision depends on which model comes next, and it is not something to fix blindly.

Summary of what I would do

To pull the three problems together:

Centre and scale everything. This is the baseline step for any method that depends on distances between points, and it costs nothing.
Convert Ba and Fe into present or absent indicators. They are mostly zeros, no power transformation helps them, and the presence of barium is close to a direct marker for Type 7.
Apply a log or square root to K, Ca and Fe. The square root brings K from 6.46 down to 0.86, and a log brings Ca from 2.02 to 1.15.
Investigate rows 172 and 173 rather than transforming them away. They look like near duplicate records and they alone drive most of the skewness in K. If they are genuine duplicates the right fix is to remove one, not to reshape the whole column.
Handle the RI and Ca correlation based on the model chosen. Drop one or use principal components for a model that is sensitive to correlation, and leave both in for a tree based model.

The general lesson from this exercise is that the right transformation depends on why a predictor looks odd. K is skewed because of two extreme rows, Ba is skewed because it is mostly zeros, and Ca is skewed in the ordinary way. Those three causes look similar on a skewness table but they call for three different fixes.

Exercise 3.2

Question. 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 (for example temperature, precipitation) and plant conditions (for example leaf spots, mold growth). The outcome labels consist of 19 distinct classes. The data can be loaded via:

library(mlbench)
data(Soybean)
## See ?Soybean for details

1. Investigate the frequency distributions for the categorical predictors. Are any of the distributions degenerate in the ways discussed earlier in this chapter?
1. Roughly 18 percent of the data are missing. Are there particular predictors that are more likely to be missing? Is the pattern of missing data related to the classes?
1. Develop a strategy for handling missing data, either by eliminating predictors or imputation.

My thought process. This dataset is the opposite of the Glass data in almost every way. The predictors are categories rather than numbers, so skewness does not apply and I cannot draw histograms of values. Instead the thing that goes wrong with a categorical predictor is that one category takes over almost everything, which leaves the predictor with almost no information in it.

The missing data is the bigger issue here, and part b is really asking whether the gaps are random or structured. That distinction decides everything in part c, because random gaps can reasonably be filled in but structured gaps usually cannot.

data(Soybean)

dim(Soybean)

## [1] 683  36

nlevels(Soybean$Class)

## [1] 19

(a) Frequency distributions and degenerate predictors

My thought process. A degenerate predictor is one that is nearly constant, so it cannot help tell the classes apart. For a categorical predictor the way to spot this is to look at how lopsided the category counts are. The caret package has nearZeroVar() which measures exactly that, using two numbers: the ratio of the most common category to the second most common, and the percentage of distinct values. I use that rather than judging by eye, so the same rule is applied to all 35 predictors.

predictors <- Soybean[, -1]

nzv_metrics <- nearZeroVar(predictors, saveMetrics = TRUE)

# The predictors flagged as near zero variance.
nzv_metrics |>
  tibble::rownames_to_column("Predictor") |>
  filter(nzv | zeroVar) |>
  select(Predictor, freqRatio, percentUnique, zeroVar, nzv) |>
  mutate(freqRatio = round(freqRatio, 2))

##   Predictor freqRatio percentUnique zeroVar  nzv
## 1 leaf.mild     26.75     0.4392387   FALSE TRUE
## 2  mycelium    106.50     0.2928258   FALSE TRUE
## 3 sclerotia     31.25     0.2928258   FALSE TRUE

Are any distributions degenerate? Yes, three of them. The function flags leaf.mild, mycelium and sclerotia.

# The actual counts behind those three flags.
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

The tables show why they were flagged. For mycelium, 639 samples sit in category 0 and only 6 sit in category 1, which is a frequency ratio of about 106 to 1. For sclerotia it is 625 against 20, and for leaf.mild it is 535 against 20 and 20.

In each case one category covers well over 90 percent of the data. A predictor like that is almost constant, so for nearly every sample it says the same thing and cannot help separate the classes.

# Visual check on a set of predictors, including the three flagged ones.
show_these <- c("mycelium", "sclerotia", "leaf.mild",
                "int.discolor", "leaf.malf", "precip",
                "temp", "plant.stand", "roots")

Soybean |>
  select(all_of(show_these)) |>
  # Some predictors are ordered factors and some are plain factors, so convert
  # them all to text before stacking them into one column.
  mutate(across(everything(), as.character)) |>
  pivot_longer(everything(), names_to = "Predictor", values_to = "Level") |>
  mutate(Level = tidyr::replace_na(Level, "missing")) |>
  ggplot(aes(x = Level)) +
  geom_bar(fill = "steelblue") +
  facet_wrap(~ Predictor, scales = "free", ncol = 3) +
  labs(title = "Frequency distribution of selected predictors",
       subtitle = "Missing values are shown as their own bar",
       x = "Category level", y = "Count")

The plot makes the contrast clear. The three flagged predictors, mycelium, sclerotia and leaf.mild, are each one tall bar with almost nothing beside it. Compare those with precip, temp and plant.stand, which are lopsided too but still have real counts spread across more than one category, so they carry usable information.

I plotted the missing values as their own bar rather than dropping them, because for this dataset the missingness turns out to be a major finding in its own right. That is what part b looks at next.

One caution about dropping them. A near zero variance flag is a warning, not an instruction. A rare category can still be very informative if it lines up with a rare class. Since some classes here have only 8 or 14 samples, a predictor that is positive in just 6 cases could in principle pick out one of those small classes. I would check that before discarding these three rather than removing them automatically.

(b) Which predictors are missing, and is it related to class?

My thought process. The question already tells me a lot is missing, so counting the total is not the interesting part. What matters is the pattern. I look at it from two directions, first by predictor and then by class, because those answer different questions and the second one turned out to be far more revealing.

# Two different ways to state how much is missing.
c(percent_of_cells_missing = round(mean(is.na(Soybean[, -1])) * 100, 2),
  percent_of_rows_affected = round(mean(!complete.cases(Soybean)) * 100, 2))

## percent_of_cells_missing percent_of_rows_affected 
##                     9.78                    17.72

A note on the 18 percent. The question says roughly 18 percent of the data are missing. That figure refers to rows, not to individual values. About 17.7 percent of the rows have at least one missing entry, while only about 9.8 percent of the individual predictor values are actually missing. Both numbers are correct, they just measure different things, and it is worth being clear about which one is meant.

missing_by_predictor <- predictors |>
  summarise(across(everything(), \(x) sum(is.na(x)))) |>
  pivot_longer(everything(), names_to = "Predictor", values_to = "Missing") |>
  mutate(Percent = round(100 * Missing / nrow(Soybean), 1)) |>
  arrange(desc(Missing))

head(missing_by_predictor, 15)

## # A tibble: 15 × 3
##    Predictor       Missing Percent
##    <chr>             <int>   <dbl>
##  1 hail                121    17.7
##  2 sever               121    17.7
##  3 seed.tmt            121    17.7
##  4 lodging             121    17.7
##  5 germ                112    16.4
##  6 leaf.mild           108    15.8
##  7 fruiting.bodies     106    15.5
##  8 fruit.spots         106    15.5
##  9 seed.discolor       106    15.5
## 10 shriveling          106    15.5
## 11 leaf.shread         100    14.6
## 12 seed                 92    13.5
## 13 mold.growth          92    13.5
## 14 seed.size            92    13.5
## 15 leaf.halo            84    12.3

Are particular predictors more likely to be missing? Yes. Four predictors are missing 121 times each, which is 17.7 percent of all rows: hail, sever, seed.tmt and lodging, with germ just behind at 112. After those comes a group including leaf.mild, fruiting.bodies, fruit.spots, seed.discolor and shriveling, each missing over 100 times.

# How many predictors have no missing values at all?
sum(missing_by_predictor$Missing == 0)

## [1] 1

missing_by_predictor |> filter(Missing == 0)

## # A tibble: 1 × 3
##   Predictor Missing Percent
##   <chr>       <int>   <dbl>
## 1 leaves          0       0

Only one predictor, leaves, is complete. Every other predictor has at least some missing values.

Is the pattern related to the classes?

missing_by_class <- Soybean |>
  group_by(Class) |>
  summarise(
    n_samples     = n(),
    missing_cells = sum(is.na(across(everything()))),
    .groups = "drop"
  ) |>
  mutate(percent_missing = round(100 * missing_cells / (n_samples * 35), 1)) |>
  arrange(desc(percent_missing))

missing_by_class

## # A tibble: 19 × 4
##    Class                       n_samples missing_cells percent_missing
##    <fct>                           <int>         <int>           <dbl>
##  1 2-4-d-injury                       16           450            80.4
##  2 cyst-nematode                      14           336            68.6
##  3 herbicide-injury                    8           160            57.1
##  4 phytophthora-rot                   88          1214            39.4
##  5 diaporthe-pod-&-stem-blight        15           177            33.7
##  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             0

Yes, and this is the most important result in the whole exercise. The missing data is not spread evenly across the classes. It is concentrated entirely in five of the nineteen classes, and the other fourteen classes have no missing values at all.

The five affected classes are 2-4-d-injury at about 80 percent missing, cyst-nematode at 69 percent, herbicide-injury at 57 percent, phytophthora-rot at 39 percent, and diaporthe-pod-&-stem-blight at 34 percent.

# Within the worst classes, is every row affected?
worst_three <- c("2-4-d-injury", "cyst-nematode", "herbicide-injury")

Soybean |>
  filter(Class %in% worst_three) |>
  mutate(na_in_row = rowSums(is.na(across(-Class)))) |>
  group_by(Class) |>
  summarise(
    n_samples           = n(),
    rows_with_any_na    = sum(na_in_row > 0),
    avg_missing_per_row = round(mean(na_in_row), 1),
    .groups = "drop"
  )

## # A tibble: 3 × 4
##   Class            n_samples rows_with_any_na avg_missing_per_row
##   <fct>                <int>            <int>               <dbl>
## 1 2-4-d-injury            16               16                28.1
## 2 cyst-nematode           14               14                24  
## 3 herbicide-injury         8                8                20

This confirms the gaps are structural rather than accidental. In 2-4-d-injury all 16 rows are affected, and each row is missing about 28 of its 35 predictors. It is not that a few values were lost here and there. Almost the entire record is absent for every sample in that class.

# For 2-4-d-injury, which predictors were actually recorded?
sub_24d <- Soybean |> filter(Class == "2-4-d-injury") |> select(-Class)

names(sub_24d)[colSums(is.na(sub_24d)) == 0]

## [1] "leaves"    "leaf.halo" "leaf.marg" "leaf.size" "leaf.malf"

This is the clearest evidence of all. For 2-4-d-injury the only predictors recorded are leaves, leaf.halo, leaf.marg, leaf.size and leaf.malf. Every one of those is about the leaves.

That tells me the missingness is not random at all. It looks like these samples were diagnosed by looking at the leaves alone, so the other measurements were never taken in the first place. Two of the affected classes, 2-4-d-injury and herbicide-injury, are chemical damage rather than disease, which fits the idea. If the damage can be identified from the leaves, there is no reason to examine the roots, stem, seed and fruit.

Why this matters so much. There is a real difference between a value that exists but was not written down and a value that was never measured because it was not relevant. The first can be estimated from similar samples. The second cannot, because there is nothing to estimate. It also means the missingness itself carries information, since knowing that a sample has most fields blank already tells you it is likely one of those five classes.

(c) A strategy for handling the missing data

My thought process. Now that I know the gaps are structural rather than random, the usual advice to impute does not fit cleanly. Imputing a value for a measurement that was never taken means inventing data. So rather than picking one method up front, I test the obvious options, report what each one actually achieves, and then give the recommendation I would follow.

# Option 1: drop the predictors with the most missing values.
drop_worst_five <- c("hail", "sever", "seed.tmt", "lodging", "germ")

after_dropping <- Soybean |> select(-all_of(drop_worst_five))

c(complete_rows_before = sum(complete.cases(Soybean)),
  complete_rows_after  = sum(complete.cases(after_dropping)),
  total_rows           = nrow(Soybean))

## complete_rows_before  complete_rows_after           total_rows 
##                  562                  562                  683

Option 1: drop the worst predictors. This is the option that looks most obvious and it turns out to achieve very little. Removing the five predictors with the most missing values barely changes the number of complete rows.

The reason is that the same rows are missing many predictors at once. Taking away five columns does not rescue a row that is still missing twenty others. This is exactly why it is worth testing an idea rather than assuming it works.

# Option 2: set aside the five classes where the missingness is concentrated.
affected_classes <- c("2-4-d-injury", "cyst-nematode", "herbicide-injury",
                      "phytophthora-rot", "diaporthe-pod-&-stem-blight")

without_affected <- Soybean |> filter(!Class %in% affected_classes)

c(rows_kept     = nrow(without_affected),
  complete_rows = sum(complete.cases(without_affected)),
  rows_lost     = nrow(Soybean) - nrow(without_affected))

##     rows_kept complete_rows     rows_lost 
##           542           542           141

Option 2: separate out the affected classes. Removing those five classes leaves 542 rows with no missing values whatsoever. The remaining data is perfectly complete, which confirms again that the missingness follows the class rather than being scattered at random.

The cost is real though. It throws away 141 samples and five of the nineteen diseases, so the model could never diagnose those conditions at all. For a tool meant to identify soybean disease, quietly dropping five diseases would be a serious limitation.

Option 3: impute. For the fourteen clean classes there is nothing to impute. For the five affected ones, imputation would mean guessing at measurements that were never taken. A method such as k nearest neighbours fills a gap using whichever samples look most similar, but with 28 of 35 fields absent there is very little left to match on, so the filled in values would mostly reflect the method rather than the plant.

Option 4: treat missing as its own category. Since these predictors are already categorical, missing can simply become another level rather than a gap that has to be filled.

What I would actually do

My recommendation is option 4, with parts of the others alongside it.

Add “not recorded” as an explicit category for each categorical predictor. This is my main recommendation. It invents no values, it keeps all 683 samples and all 19 classes, and it turns the missingness into a usable signal. Since a blank record is strongly associated with those five classes, a model can learn from the absence itself. The pattern of what was not measured is genuinely informative here.
Use a model that handles missing values natively. Tree based methods such as random forests or CART can split on a missing category directly, which pairs well with the step above. That suits this data better than a method that demands a complete matrix.
Reconsider the three degenerate predictors from part a, but check them against the small classes first rather than dropping them automatically.
Keep the affected classes in the data. Removing them gives a cleaner table and a worse tool. Those five conditions are exactly the ones a diagnostic model needs to handle, and the imbalance in class sizes is already a bigger problem than the missing values are.
Report the limitation honestly. Any model will be weaker on the five affected classes, because there is genuinely less information about them. That is a property of how the data was collected, not something a cleverer method can repair.

The general lesson from this exercise is that the right way to handle missing data depends on why it is missing. Here the gaps are not accidents, they are a record of which examinations were carried out. Filling them in would erase that information. Marking them explicitly keeps it.

Homework 4: Applied Predictive Modeling Exercises 3.1 and 3.2

Umais Siddiqui

2026-10-02