Predictive modeling

Exercises

3.1

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

#load the library
library(mlbench) 

Warning: package 'mlbench' was built under R version 4.5.2

library(tidyverse)

Warning: package 'ggplot2' was built under R version 4.5.2

Warning: package 'tibble' was built under R version 4.5.2

Warning: package 'tidyr' was built under R version 4.5.2

Warning: package 'readr' was built under R version 4.5.2

Warning: package 'purrr' was built under R version 4.5.2

Warning: package 'dplyr' was built under R version 4.5.2

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.0     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.2     ✔ tibble    3.3.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggplot2)
#data
data(Glass) 
#String
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 ...

# Summary of all variables
summary(Glass)

       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          Type  
 Min.   :0.00000   1:70  
 1st Qu.:0.00000   2:76  
 Median :0.00000   3:17  
 Mean   :0.05701   5:13  
 3rd Qu.:0.10000   6: 9  
 Max.   :0.51000   7:29  

###Exploring distribution of all predictors

predictors <- Glass[, names(Glass) != "Type"]

par(mfrow = c(3, 3))

for (name in names(predictors)) {
  hist(predictors[[name]],
       main = paste("Distribution of", name),
       xlab = name)
}

##Interpretation:

Looking at the distributions of the predictors, RI, K, Ca, Ba, and Fe appear to be right-skewed. Mg appears to be left-skewed, with most observations concentrated at the higher end of its range. Si also shows some left skewness. Na and Al appear more symmetric compared with the other predictors, although Al may have a slight right skew. Ba and Fe have many observations near zero and show particularly strong right skewness.

###Relationship between predictors

pairs(predictors)

round(cor(predictors), 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

##Interpretation:

The scatterplot matrix and correlation matrix show relationships among several predictors. The strongest relationship is between RI and Ca, with a correlation of 0.81, indicating a strong positive relationship. RI and Si have a moderate negative relationship with a correlation of -0.54. Other moderate relationships include Mg and Ba (-0.49), Mg and Al (-0.48), Al and Ba (0.48), and Mg and Ca (-0.44). Many of the remaining predictor pairs have relatively weak linear relationships.

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

par(mfrow = c(3, 3))

for (name in names(predictors)) {
  boxplot(predictors[[name]],
          main = name)
}

##Interpretation:

Based on the boxplots, most predictors appear to contain potential outliers. RI, Na, Al, Si, K, Ca, Ba, and Fe have observations outside the boxplot whiskers, while Mg does not show obvious outliers. From the histograms, RI, K, Ca, Ba, and Fe are right-skewed, while Mg and Si show left skewness. Na and Al appear more symmetric compared with the other predictors.

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

par(mfrow = c(2, 2))

hist(Glass$K, main = "Original K")
hist(sqrt(Glass$K), main = "Square Root")
hist(log1p(Glass$K), main = "Log(1 + K)")
hist(log(Glass$K + 0.01), main = "Log(K + 0.01)")

par(mfrow = c(2, 2))

hist(Glass$Ca, main = "Original Ca")
hist(sqrt(Glass$Ca), main = "Square Root")
hist(log1p(Glass$Ca), main = "Log(1 + Ca)")
hist(log(Glass$Ca), main = "Log(Ca)")

par(mfrow = c(2, 2))

hist(Glass$Ba, main = "Original Ba")
hist(sqrt(Glass$Ba), main = "Square Root")
hist(log1p(Glass$Ba), main = "Log(1 + Ba)")
hist(log(Glass$Ba + 0.01), main = "Log(Ba + 0.01)")

par(mfrow = c(2, 2))

hist(Glass$Fe, main = "Original Fe")
hist(sqrt(Glass$Fe), main = "Square Root")
hist(log1p(Glass$Fe), main = "Log(1 + Fe)")
hist(log(Glass$Fe + 0.01), main = "Log(Fe + 0.01)")

##Interpretation:

Several predictors are skewed, so transformations may be considered. For Ca, square-root and log transformations reduce the long right tail, with a log transformation providing a more compressed distribution. K is also strongly right-skewed, but it contains many zero values, making a standard log transformation inappropriate. Square-root and log1p transformations can accommodate the zeros but do not completely remove the skewness. Ba and Fe contain a large number of zero values, so their skewness is partly caused by the concentration of observations at zero. Log transformations with an added constant preserve these observations but do not necessarily produce symmetric distributions. Therefore, transformations may be useful for some predictors, especially Ca, but their effect on classification performance should ultimately be evaluated during model training.

3.2

data(Soybean)
#Summary
summary(Soybean)

                 Class          date     plant.stand  precip      temp    
 brown-spot         : 92   5      :149   0   :354    0   : 74   0   : 80  
 alternarialeaf-spot: 91   4      :131   1   :293    1   :112   1   :374  
 frog-eye-leaf-spot : 91   3      :118   NA's: 36    2   :459   2   :199  
 phytophthora-rot   : 88   2      : 93               NA's: 38   NA's: 30  
 anthracnose        : 44   6      : 90                                    
 brown-stem-rot     : 44   (Other):101                                    
 (Other)            :233   NA's   :  1                                    
   hail     crop.hist  area.dam    sever     seed.tmt     germ     plant.growth
 0   :435   0   : 65   0   :123   0   :195   0   :305   0   :165   0   :441    
 1   :127   1   :165   1   :227   1   :322   1   :222   1   :213   1   :226    
 NA's:121   2   :219   2   :145   2   : 45   2   : 35   2   :193   NA's: 16    
            3   :218   3   :187   NA's:121   NA's:121   NA's:112               
            NA's: 16   NA's:  1                                                
                                                                               
                                                                               
 leaves  leaf.halo  leaf.marg  leaf.size  leaf.shread leaf.malf  leaf.mild 
 0: 77   0   :221   0   :357   0   : 51   0   :487    0   :554   0   :535  
 1:606   1   : 36   1   : 21   1   :327   1   : 96    1   : 45   1   : 20  
         2   :342   2   :221   2   :221   NA's:100    NA's: 84   2   : 20  
         NA's: 84   NA's: 84   NA's: 84                          NA's:108  
                                                                           
                                                                           
                                                                           
   stem     lodging    stem.cankers canker.lesion fruiting.bodies ext.decay 
 0   :296   0   :520   0   :379     0   :320      0   :473        0   :497  
 1   :371   1   : 42   1   : 39     1   : 83      1   :104        1   :135  
 NA's: 16   NA's:121   2   : 36     2   :177      NA's:106        2   : 13  
                       3   :191     3   : 65                      NA's: 38  
                       NA's: 38     NA's: 38                                
                                                                            
                                                                            
 mycelium   int.discolor sclerotia  fruit.pods fruit.spots   seed    
 0   :639   0   :581     0   :625   0   :407   0   :345    0   :476  
 1   :  6   1   : 44     1   : 20   1   :130   1   : 75    1   :115  
 NA's: 38   2   : 20     NA's: 38   2   : 14   2   : 57    NA's: 92  
            NA's: 38                3   : 48   4   :100              
                                    NA's: 84   NA's:106              
                                                                     
                                                                     
 mold.growth seed.discolor seed.size  shriveling  roots    
 0   :524    0   :513      0   :532   0   :539   0   :551  
 1   : 67    1   : 64      1   : 59   1   : 38   1   : 86  
 NA's: 92    NA's:106      NA's: 92   NA's:106   2   : 15  
                                                 NA's: 31  
                                                           
                                                           
                                                           

#Soybean

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

#str(Soybean)
#summary(Soybean)
Soybean_long <- Soybean |>
  mutate(across(-Class, as.character)) |>
  pivot_longer(
    cols = -Class,
    names_to = "Predictor",
    values_to = "Value"
  )
Soybean_long

# A tibble: 23,905 × 3
   Class                 Predictor   Value
   <fct>                 <chr>       <chr>
 1 diaporthe-stem-canker date        6    
 2 diaporthe-stem-canker plant.stand 0    
 3 diaporthe-stem-canker precip      2    
 4 diaporthe-stem-canker temp        1    
 5 diaporthe-stem-canker hail        0    
 6 diaporthe-stem-canker crop.hist   1    
 7 diaporthe-stem-canker area.dam    1    
 8 diaporthe-stem-canker sever       1    
 9 diaporthe-stem-canker seed.tmt    0    
10 diaporthe-stem-canker germ        0    
# ℹ 23,895 more rows

Soybean_frequency <- Soybean_long |>
  count(Predictor, Value)
Soybean_frequency

# A tibble: 133 × 3
   Predictor     Value     n
   <chr>         <chr> <int>
 1 area.dam      0       123
 2 area.dam      1       227
 3 area.dam      2       145
 4 area.dam      3       187
 5 area.dam      <NA>      1
 6 canker.lesion 0       320
 7 canker.lesion 1        83
 8 canker.lesion 2       177
 9 canker.lesion 3        65
10 canker.lesion <NA>     38
# ℹ 123 more rows

Soybean_long |>
  filter(!is.na(Value)) |>
  ggplot(aes(x = Value)) +
  geom_bar() +
  facet_wrap(~ Predictor, scales = "free_x") +
  labs(
    title = "Frequency Distribution of Soybean Predictor",
    x = "Category",
    y = "Frequency"
  ) +
  theme_minimal()

Soybean_long |>
  filter(!is.na(Value)) |>
  count(Predictor, Value) |>
  group_by(Predictor) |>
  mutate(Percent = n / sum(n) * 100) |>
  summarise(
    Largest_Percent = max(Percent),
    Number_of_Categories = n()
  ) |>
  arrange(desc(Largest_Percent))

# A tibble: 35 × 3
   Predictor     Largest_Percent Number_of_Categories
   <chr>                   <dbl>                <int>
 1 mycelium                 99.1                    2
 2 sclerotia                96.9                    2
 3 shriveling               93.4                    2
 4 leaf.mild                93.0                    3
 5 lodging                  92.5                    2
 6 leaf.malf                92.5                    2
 7 int.discolor             90.1                    3
 8 seed.size                90.0                    2
 9 seed.discolor            88.9                    2
10 leaves                   88.7                    2
# ℹ 25 more rows

##Interpretation:

The frequency distributions show that several categorical predictors are highly unbalanced. For example, mycelium is dominated by category 0, with 639 observations compared with only 6 observations in category 1. Other predictors such as sclerotia, leaf.malf, shriveling, and leaf.mild also show strong imbalance toward one category. In contrast, variables such as plant.stand and crop.hist have more balanced distributions across their categories. Therefore, some predictors appear to have near-zero variation and may provide limited information for prediction.

b. Roughly 18% 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?

Soybean_long |>
  group_by(Predictor) |>
  summarize(Missing = sum(is.na(Value)),
            Percent_Missing = mean(is.na(Value))* 100) |>
  arrange(desc(Percent_Missing))

# A tibble: 35 × 3
   Predictor       Missing Percent_Missing
   <chr>             <int>           <dbl>
 1 hail                121            17.7
 2 lodging             121            17.7
 3 seed.tmt            121            17.7
 4 sever               121            17.7
 5 germ                112            16.4
 6 leaf.mild           108            15.8
 7 fruit.spots         106            15.5
 8 fruiting.bodies     106            15.5
 9 seed.discolor       106            15.5
10 shriveling          106            15.5
# ℹ 25 more rows

##Let examine which predictors are missing within each class

Soybean_long |>
  group_by(Class, Predictor) |>
  summarise(Percent_Missing = mean(is.na(Value)) * 100,
            .groups = "drop") |>
  filter(Percent_Missing > 0) |>
  arrange(Class, desc(Percent_Missing))

# A tibble: 106 × 3
   Class        Predictor       Percent_Missing
   <fct>        <chr>                     <dbl>
 1 2-4-d-injury canker.lesion               100
 2 2-4-d-injury crop.hist                   100
 3 2-4-d-injury ext.decay                   100
 4 2-4-d-injury fruit.pods                  100
 5 2-4-d-injury fruit.spots                 100
 6 2-4-d-injury fruiting.bodies             100
 7 2-4-d-injury germ                        100
 8 2-4-d-injury hail                        100
 9 2-4-d-injury int.discolor                100
10 2-4-d-injury leaf.mild                   100
# ℹ 96 more rows

##Interpretation:

There are particular predictors that are more likely to have missing values than others. The missing values are also not evenly distributed across the classes. Some classes, such as 2-4-d-injury, cyst-nematode, and herbicide-injury, have several predictors with 100% missing values, while phytophthora-rot has many predictors with about 77% missing values. Therefore, the pattern of missing data appears to be related to the classes.

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

get_mode <- function(x) {
  values <- x[!is.na(x)]
  names(sort(table(values), decreasing = TRUE))[1]
}

Soybean_imputed <- Soybean |>
  mutate(across(-Class,
                ~ replace(., is.na(.), get_mode(.))
                ))
#
head(Soybean_imputed, 10)

                   Class date plant.stand precip temp hail crop.hist area.dam
1  diaporthe-stem-canker    6           0      2    1    0         1        1
2  diaporthe-stem-canker    4           0      2    1    0         2        0
3  diaporthe-stem-canker    3           0      2    1    0         1        0
4  diaporthe-stem-canker    3           0      2    1    0         1        0
5  diaporthe-stem-canker    6           0      2    1    0         2        0
6  diaporthe-stem-canker    5           0      2    1    0         3        0
7  diaporthe-stem-canker    5           0      2    1    0         2        0
8  diaporthe-stem-canker    4           0      2    1    1         1        0
9  diaporthe-stem-canker    6           0      2    1    0         3        0
10 diaporthe-stem-canker    4           0      2    1    0         2        0
   sever seed.tmt germ plant.growth leaves leaf.halo leaf.marg leaf.size
1      1        0    0            1      1         0         2         2
2      2        1    1            1      1         0         2         2
3      2        1    2            1      1         0         2         2
4      2        0    1            1      1         0         2         2
5      1        0    2            1      1         0         2         2
6      1        0    1            1      1         0         2         2
7      1        1    0            1      1         0         2         2
8      1        0    2            1      1         0         2         2
9      1        1    1            1      1         0         2         2
10     2        0    2            1      1         0         2         2
   leaf.shread leaf.malf leaf.mild stem lodging stem.cankers canker.lesion
1            0         0         0    1       1            3             1
2            0         0         0    1       0            3             1
3            0         0         0    1       0            3             0
4            0         0         0    1       0            3             0
5            0         0         0    1       0            3             1
6            0         0         0    1       0            3             0
7            0         0         0    1       1            3             1
8            0         0         0    1       0            3             1
9            0         0         0    1       0            3             1
10           0         0         0    1       0            3             1
   fruiting.bodies ext.decay mycelium int.discolor sclerotia fruit.pods
1                1         1        0            0         0          0
2                1         1        0            0         0          0
3                1         1        0            0         0          0
4                1         1        0            0         0          0
5                1         1        0            0         0          0
6                1         1        0            0         0          0
7                1         1        0            0         0          0
8                1         1        0            0         0          0
9                1         1        0            0         0          0
10               1         1        0            0         0          0
   fruit.spots seed mold.growth seed.discolor seed.size shriveling roots
1            4    0           0             0         0          0     0
2            4    0           0             0         0          0     0
3            4    0           0             0         0          0     0
4            4    0           0             0         0          0     0
5            4    0           0             0         0          0     0
6            4    0           0             0         0          0     0
7            4    0           0             0         0          0     0
8            4    0           0             0         0          0     0
9            4    0           0             0         0          0     0
10           4    0           0             0         0          0     0

sum(is.na(Soybean_imputed))

[1] 0

##Interpretation:

I would use imputation rather than eliminate predictors because removing predictors could result in the loss of useful information. Since the predictors are categorical, I used mode imputation, replacing each missing value with the most frequent observed category for that predictor. After imputation, there were no remaining missing values. However, because the missing-data pattern appears to be related to the disease classes, mode imputation should be used carefully, especially for predictors with high percentages of missing values.