Let’s first load the library and open the data. I will install the mlbench package in my console so that the system won’t have to do it over every time that i run my entire file.
library(mlbench)
Warning: package 'mlbench' was built under R version 4.5.3
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 ...
Our Glass data have nine predictors including RI(refractive index), Na,Mg,Al,Si,K, CA, Ba and Fe.
Using visualizations,let’s explore the predictor variables to understand their distributions as well as the relationship between predictors.
One of the best way to explore variables’ distributions is through an histogram. In order to create the histogram for our predictor, we will have to extract only the predictors from our Glass data meaning we should remove the variable “type” from the data.
# Let's create our dataset of interestpredictors <- Glass[, -10] # drop Typehead(predictors)
RI Na Mg Al Si K Ca Ba Fe
1 1.52101 13.64 4.49 1.10 71.78 0.06 8.75 0 0.00
2 1.51761 13.89 3.60 1.36 72.73 0.48 7.83 0 0.00
3 1.51618 13.53 3.55 1.54 72.99 0.39 7.78 0 0.00
4 1.51766 13.21 3.69 1.29 72.61 0.57 8.22 0 0.00
5 1.51742 13.27 3.62 1.24 73.08 0.55 8.07 0 0.00
6 1.51596 12.79 3.61 1.62 72.97 0.64 8.07 0 0.26
# Let's create a histogram for each predictorfor (var innames(predictors)) hist(predictors[[var]], main = var, xlab ="")
# Let's find relationships between predictors by looking at the correlation among predictorslibrary(corrplot) # Loading the coorplot library after installing it in my console
Warning: package 'corrplot' was built under R version 4.5.3
corrplot 0.95 loaded
# Let's create a correlation matrix and the plot associatedcor_matrix <-round(cor(predictors),3) # generating the correlation matrix at 3 decimals placescor_matrix
RI Na Mg Al Si K Ca Ba Fe
RI 1.000 -0.192 -0.122 -0.407 -0.542 -0.290 0.810 0.000 0.143
Na -0.192 1.000 -0.274 0.157 -0.070 -0.266 -0.275 0.327 -0.241
Mg -0.122 -0.274 1.000 -0.482 -0.166 0.005 -0.444 -0.492 0.083
Al -0.407 0.157 -0.482 1.000 -0.006 0.326 -0.260 0.479 -0.074
Si -0.542 -0.070 -0.166 -0.006 1.000 -0.193 -0.209 -0.102 -0.094
K -0.290 -0.266 0.005 0.326 -0.193 1.000 -0.318 -0.043 -0.008
Ca 0.810 -0.275 -0.444 -0.260 -0.209 -0.318 1.000 -0.113 0.125
Ba 0.000 0.327 -0.492 0.479 -0.102 -0.043 -0.113 1.000 -0.059
Fe 0.143 -0.241 0.083 -0.074 -0.094 -0.008 0.125 -0.059 1.000
corrplot(cor_matrix, order ="hclust") # we are plotting the correlation using hierarchical clustering order.
# Let's create a pairwise scatterplot matrix to see the relationship between predictorspairs(predictors)
Interpretation
When we look at the correlation matrix, the correlation coefficient between RI and Ca is equal to 0.810 which means there is a strong positive relationship between these two predictors. We also observe some moderate relationship between all predictors with correlation coefficient between (+/-) 0.30 and 0.69 and finally we have some coefficients less than (+/-) 0.3 which shows some weak relationship between some predictors.
Let’s find any possible outliers or skewness in predictors.
The box plot could easily help us figure out any outliers. In regard to the skewness in predictors,
# Boxplots flag outliersfor (var innames(predictors)) boxplot(predictors[[var]], main = var)
Interpretation:
Following all the box plot, we observed outliers in predictors K, Ba, Fe, Ca and Al with some observed high values in the respective graph with some value that are single out and not directly related to the majority of the data.
When visualizing all the histograms form the question a), overall the Na and Al predictors histogram are roughly symmetrical and exhibit a sort of bell-shape form. As for the K, Ba, Ca, Fe, and Ri predictor exhibit a long tail to the right which mean they are right skewed. The Mg predictor shows what can be considered like 2 peaks which makes it a bimodal distribution. Lastly, the Si predictor is mainly symmetrical but also exhibit a long tail to the left.
Are there any relevant transformations of one or more predictors that might improved the classification model?
Since we discover some skew distribution, some outliers and zeros in the predictors, therefore we believe some transformations can be done to improve the classification model. The necessary transformations are as followed:
We can use the Yeo-Johnson transformation handles zero and negative values and for the skewed predictors.
We can use the Spatial sign transformation to deal with outliers.
We can use the principal component analysis to deal with correlated predictors.
# We install and load the caret to deal with the pre-processing steps in R. I installed the package in my console to avoid for it to be installed anytime we run our command. library(caret)
Warning: package 'caret' was built under R version 4.5.3
Loading required package: ggplot2
Warning: package 'ggplot2' was built under R version 4.5.3
Loading required package: lattice
Warning: package 'lattice' was built under R version 4.5.2
# Box-Cox needs strictly positive values, so it works for RI, Na, Al, Si, Ca# Yeo-Johnson transformation for (Mg, K, Ba, Fe ) that contain zerospreProcess_Data <-preProcess(predictors, method =c("YeoJohnson", "center", "scale"))YeoJohnson_trans <-predict(preProcess_Data, predictors)# Spatial sign transformation to shrinks outlierspreProcess2_Data <-preProcess(predictors, method =c("center", "scale"))ss_trans <-spatialSign(predict(preProcess2_Data, predictors))head(ss_trans)
RI Na Mg Al Si K
1 0.3862138 0.12608198 0.5551342 -0.3063816 -0.49869486 -0.29720647
2 -0.1782866 0.42318047 0.4548937 -0.1218880 0.07316353 -0.01874385
3 -0.4739625 0.09851782 0.3951817 0.1254439 0.28831726 -0.10811108
4 -0.1932662 -0.20158533 0.5799799 -0.2581474 -0.04397201 0.09305653
5 -0.2268925 -0.12303180 0.4726731 -0.2991170 0.40373539 0.05916450
6 -0.3036441 -0.29004891 0.2459643 0.1344417 0.15791804 0.08402131
Ca Ba Fe
1 -0.06449657 -0.1561358 -0.2594843
2 -0.56756147 -0.2523255 -0.4193433
3 -0.54468469 -0.2318677 -0.3853442
4 -0.43085072 -0.2929133 -0.4867968
5 -0.45422808 -0.2565822 -0.4264175
6 -0.23892010 -0.1349600 0.7986264
# Principal Component Analysis to handle correlated predictors such as RI and CapreProcess3_Data <-preProcess(predictors, method =c("YeoJohnson", "center", "scale", "pca"))preProcess3_Data
Created from 214 samples and 9 variables
Pre-processing:
- centered (9)
- ignored (0)
- principal component signal extraction (9)
- scaled (9)
- Yeo-Johnson transformation (5)
Lambda estimates for Yeo-Johnson transformation:
-0.19, 2.24, 0, -0.98, -1.36
PCA needed 6 components to capture 95 percent of the variance
Let’s investigate the frequency distribution for the categorical predictors. Are any of the distribution degenerate in the way discussed earlier in this chapter?
Based on our chapter 3, degenerate distributions refer to predictors that have zero or near-zero variance(NZV). In order to check that we can use the caret package.
library(caret)# Frequency tables for every predictorlapply(Soybean[, -1], table, useNA ="ifany")
# Let's check for zero and near-zero-variance predictorsnearzerovar_info <-nearZeroVar(predictors, saveMetrics =TRUE)nearzerovar_info[nearzerovar_info$nearzerov ==TRUE, ]
# Visualize or tabulate one of the degenerate variables (e.g., leaf.malf, mycelium)table(predictors$leaf.malf, useNA ="always")
<NA>
0
table(predictors$mycelium, useNA ="always")
<NA>
0
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?
First, we will check the percent missing by predictor.
# Let's generate the missing data by class tablemissingdata_by_class <-aggregate( . ~ Class,data =cbind(Class = Soybean$Class, as.data.frame(is.na(Soybean[, -1]))),FUN = mean)head(missingdata_by_class,10)
