DATA624

Exercise 3.1

The UC Irvine Machine Learning Repository6 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.

I couldn’t get this ‘Glass’ dataset to load without a fatal error so I had to load it from the original source ‘evclass’

glass_hw <- as.data.frame(evclass::glass)
names(glass_hw) <- c("RI", "Na", "Mg", "Al", "Si", "K", "Ca", "Ba", "Fe", "Type")

str(glass_hw)

## 'data.frame':    185 obs. of  10 variables:
##  $ RI  : num  1.52 1.52 1.52 1.52 1.52 ...
##  $ Na  : num  13.4 14.2 12.9 13.9 13.6 ...
##  $ Mg  : num  3.39 3.82 3.62 3.6 3.35 2.19 2.39 2.76 1.35 3.44 ...
##  $ Al  : num  1.28 0.47 1.57 1.36 1.23 1.66 1.56 0.83 1.63 1.45 ...
##  $ Si  : num  72.6 71.8 73 72.7 72.1 ...
##  $ K   : num  0.52 0.11 0.61 0.48 0.59 0 0 0.35 0.39 0.44 ...
##  $ Ca  : num  8.65 9.57 8.11 7.83 8.91 ...
##  $ Ba  : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Fe  : num  0 0 0 0 0 0 0 0.2 0.18 0 ...
##  $ Type: num  3 1 2 1 1 4 4 2 2 2 ...

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

The first step is to visualize the distributions of all predictors. We can see that some of the distributions are very skewed.

library(Hmisc)

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:dplyr':
## 
##     src, summarize

## The following objects are masked from 'package:base':
## 
##     format.pval, units

hist.data.frame(glass_hw)

We can look at the correlations between predictors using the ‘cor’ function

correlations <- cor(glass_hw[,1:9])
dim(correlations)

## [1] 9 9

correlations[1:4, 1:4]

##            RI          Na          Mg         Al
## RI  1.0000000 -0.16626897 -0.34970064 -0.3756344
## Na -0.1662690  1.00000000  0.06220203 -0.2155130
## Mg -0.3497006  0.06220203  1.00000000 -0.2126572
## Al -0.3756344 -0.21551297 -0.21265716  1.0000000

We can also plot the correlations which is more informative. Interesting, there are a lot of negative correlations between the variables. For example Magnesium and Calcium are highly negatively correlated, I guess glass has either one or the other. The Refractive Index is also highly positively or negatively correlated with different elements. For example, Calcium is highly correlated with the Refractive Index while Silicon is highly negatively correlated.

library(corrplot)

## corrplot 0.94 loaded

corrplot(cor(glass_hw[1:9]), order = "hclust")

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

In a way, I think this dataset is a weird choice for this problem. If I had a multi-class classification problem, my first choice of algorithms would be a tree-based method that is robust to skewed or highly correlated data. So the outcome of the all this data exploration wouldn’t be as consequential.

I can tell from the skewness that their probably outliers. Because this in a muli-class classification problem. While the classes are represented as numbers, I have no reason to believe that they are numeric so looking for a linear relationship in at scatter plots of the of the predictors and classes doesn’t make a ton of sense.

library(ggplot2)
ggplot(data = glass_hw, aes(x = RI, y = Type)) + geom_point()

When we use ‘ggpairs’ we can see the outliers in both Potassium and Barium. The distribution of Potassium, Iron and Barium are highly skewed.

GGally::ggpairs(glass_hw)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

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

I’m going to use the Box-Cox transformation to see if it helps with the skew and outliers.

library(caret)

## Loading required package: lattice

lambda <- BoxCoxTrans(y = glass_hw$Type, x = glass_hw$Fe)
lambda

## Box-Cox Transformation
## 
## 185 data points used to estimate Lambda
## 
## Input data summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   1.000   2.000   1.951   2.000   4.000 
## 
## Largest/Smallest: 4 
## Sample Skewness: 0.869 
## 
## Estimated Lambda: -0.3

I can use the ‘caret’ package method ‘preProcess’ to to calculate the Box-Cox transformation and then the predict method to apply the transformation

bc_trans <- preProcess(glass_hw, method = c("BoxCox"))
transformed <- predict(bc_trans, glass_hw)

I’m a little surprised how little the Box-Cox transformation changed the distributions.

GGally::ggpairs(transformed)

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.

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

If I was doing further work with this dataset I would recode some of the variables. For most predictors a ‘0’ indicates a non-event (e.g. seed.tmt: 0 = ‘none’) and for a few variables ‘0’ means an event (e.g. hail: 0 = ‘Yes’). There is also discrepancy between the levels and what the documentation says the levels are.

library(mlbench)
data(Soybean)

It’s hard to see in the plot but for a lot of the predictors, there are a lot of non-events (‘1’) and then lower counts for the event levels (‘2’,‘3’,…). For example for the column stem.cankers, the levels are “absent(1),below-soil(2),above-s(3),ab-sec-nde(4)”. It would make sense to bin the data as absent(0), present(1)

library(Hmisc)
soy_preds <- Soybean[,2:36]
soy_preds <- lapply(soy_preds,as.numeric)
soy_preds <- as.data.frame(soy_preds)
hist.data.frame(soy_preds)

count(soy_preds, stem.cankers)

##   stem.cankers   n
## 1            1 379
## 2            2  39
## 3            3  36
## 4            4 191
## 5           NA  38

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?

There is a pattern of missingness associated with classes. The classes with low counts (under 20 cases) tend to have high levels of missingness.

library(gtsummary)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ✔ readr     2.1.5     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()    masks stats::filter()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ purrr::lift()      masks caret::lift()
## ✖ Hmisc::src()       masks dplyr::src()
## ✖ Hmisc::summarize() masks dplyr::summarize()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Soybean %>%
  mutate(across(.cols = -Class, .fns = is.na)) %>%
  # then createa summary table tabulating missing data rates
  tbl_summary(by=Class) %>%
  modify_caption("Missing Data Counts") %>%
  as_kable()

Missing Data Counts
Characteristic	2-4-d-injury N = 16	cyst-nematode N = 14	diaporthe-pod-&-stem-blight N = 15	herbicide-injury N = 8	phytophthora-rot N = 88
date	1 (6.3%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
plant.stand	16 (100%)	14 (100%)	6 (40%)	0 (0%)	0 (0%)
precip	16 (100%)	14 (100%)	0 (0%)	8 (100%)	0 (0%)
temp	16 (100%)	14 (100%)	0 (0%)	0 (0%)	0 (0%)
hail	16 (100%)	14 (100%)	15 (100%)	8 (100%)	68 (77%)
crop.hist	16 (100%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
area.dam	1 (6.3%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
sever	16 (100%)	14 (100%)	15 (100%)	8 (100%)	68 (77%)
seed.tmt	16 (100%)	14 (100%)	15 (100%)	8 (100%)	68 (77%)
germ	16 (100%)	14 (100%)	6 (40%)	8 (100%)	68 (77%)
plant.growth	16 (100%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaves	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaf.halo	0 (0%)	14 (100%)	15 (100%)	0 (0%)	55 (63%)
leaf.marg	0 (0%)	14 (100%)	15 (100%)	0 (0%)	55 (63%)
leaf.size	0 (0%)	14 (100%)	15 (100%)	0 (0%)	55 (63%)
leaf.shread	16 (100%)	14 (100%)	15 (100%)	0 (0%)	55 (63%)
leaf.malf	0 (0%)	14 (100%)	15 (100%)	0 (0%)	55 (63%)
leaf.mild	16 (100%)	14 (100%)	15 (100%)	8 (100%)	55 (63%)
stem	16 (100%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
lodging	16 (100%)	14 (100%)	15 (100%)	8 (100%)	68 (77%)
stem.cankers	16 (100%)	14 (100%)	0 (0%)	8 (100%)	0 (0%)
canker.lesion	16 (100%)	14 (100%)	0 (0%)	8 (100%)	0 (0%)
fruiting.bodies	16 (100%)	14 (100%)	0 (0%)	8 (100%)	68 (77%)
ext.decay	16 (100%)	14 (100%)	0 (0%)	8 (100%)	0 (0%)
mycelium	16 (100%)	14 (100%)	0 (0%)	8 (100%)	0 (0%)
int.discolor	16 (100%)	14 (100%)	0 (0%)	8 (100%)	0 (0%)
sclerotia	16 (100%)	14 (100%)	0 (0%)	8 (100%)	0 (0%)
fruit.pods	16 (100%)	0 (0%)	0 (0%)	0 (0%)	68 (77%)
fruit.spots	16 (100%)	14 (100%)	0 (0%)	8 (100%)	68 (77%)
seed	16 (100%)	0 (0%)	0 (0%)	8 (100%)	68 (77%)
mold.growth	16 (100%)	0 (0%)	0 (0%)	8 (100%)	68 (77%)
seed.discolor	16 (100%)	14 (100%)	0 (0%)	8 (100%)	68 (77%)
seed.size	16 (100%)	0 (0%)	0 (0%)	8 (100%)	68 (77%)
shriveling	16 (100%)	14 (100%)	0 (0%)	8 (100%)	68 (77%)
roots	16 (100%)	0 (0%)	15 (100%)	0 (0%)	0 (0%)

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

I’m going to handle missingness in two ways. I’m going to remove the classes with less than 20 events and then I’m going to replace the missing values with the median of the column. Replacing missing values with the mean is problematic because taking the mean of factor levels that aren’t truely numeric doesn’t make sense.

I tried taking the mode but couldn’t get the code to work

Soybean_cleaned <- lapply(Soybean,as.numeric)
Soybean_cleaned <- as.data.frame(Soybean_cleaned)
Soybean_cleaned <- Soybean_cleaned %>% group_by(Class) %>% filter(n() > 19) %>% ungroup()

Soybean_cleaned <- Soybean_cleaned %>% mutate(across(everything(), ~replace_na(.x, median(na.omit(.)))))

Soybean_cleaned %>%
  mutate(across(.cols = -Class, .fns = is.na)) %>%
  # then createa summary table tabulating missing data rates
  tbl_summary(by=Class) %>%
  modify_caption("Missing Data Counts") %>%
  as_kable()

Missing Data Counts
Characteristic	2 N = 91	3 N = 44	4 N = 20	5 N = 20	6 N = 92	7 N = 44	8 N = 20	11 N = 20	12 N = 20	13 N = 91	15 N = 20	16 N = 88	17 N = 20	18 N = 20	19 N = 20
date	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
plant.stand	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
precip	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
temp	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
hail	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
crop.hist	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
area.dam	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
sever	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
seed.tmt	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
germ	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
plant.growth	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaves	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaf.halo	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaf.marg	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaf.size	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaf.shread	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaf.malf	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
leaf.mild	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
stem	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
lodging	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
stem.cankers	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
canker.lesion	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
fruiting.bodies	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
ext.decay	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
mycelium	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
int.discolor	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
sclerotia	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
fruit.pods	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
fruit.spots	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
seed	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
mold.growth	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
seed.discolor	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
seed.size	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
shriveling	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
roots	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)

DATA624_HW4

William Aiken

9/29/2024

Exercise 3.1

Exercise 3.2.