library(ggplot2)
library(tidyverse)
library(corrplot)
library(caret)
library(e1071)data_624_hw04
Exercise 3.1
library(mlbench)
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 ...a.
The data is explored below.
Histograms and skew
Glass |>
select(-Type) |>
pivot_longer(cols = everything()) |>
ggplot(aes(x = value)) +
geom_histogram(bins = 30, fill = 'skyblue', color = 'black') +
facet_wrap(~ name, scales = 'free')
Inspecting the histograms for each feature shows a variety of distributions.
Aluminum (Al) approaches normal with a bit of a right skew
Barium (Ba) is very right skewed, the majority of observations have a 0 value or do not contain Barium
Calcium (Ca) approaches normal but with a bit of a right skew
Iron (Fe) is very right skewed, the majority of observations have a 0 value or do not contain Iron.
Potassium (K) is right skewed with the majority of observations having less than 1% Potassium.
Magnesium (Mg) is bimodal with peaks at 0 and between 3% and 4%
Sodium (Na) approaches normal
Refractive Index (RI) is right skewed
Silicon (Si) approaches normal
skew_values <- Glass |>
select(-Type) |>
apply(2, skewness)
skew_values |> as.data.frame() skew_values
RI 1.6027151
Na 0.4478343
Mg -1.1364523
Al 0.8946104
Si -0.7202392
K 6.4600889
Ca 2.0184463
Ba 3.3686800
Fe 1.7298107Transformations for Skew
A possible transformation for skew is the Box-Cox tranformation.
glass_boxcox_check <- Glass |>
select(-Type) |>
lapply(BoxCoxTrans)
sapply(glass_boxcox_check, function(x) x$lambda) RI Na Mg Al Si K Ca Ba Fe
-2.0 -0.1 NA 0.5 2.0 NA -1.1 NA NA The Box Cox transformation is not available for use with Mg, K, Ba, Fe. Box Cox requires strictly positive values, and those features contain 0 values. Those elements are largely the ones with the most significant skew.
glass_boxcox <- Glass
glass_boxcox$RI <- glass_boxcox$RI |> BoxCoxTrans() |> predict(newdata = Glass$RI)
glass_boxcox$Na <- glass_boxcox$Na |> BoxCoxTrans() |> predict(newdata = Glass$Na)
glass_boxcox$Al <- glass_boxcox$Al |> BoxCoxTrans() |> predict(newdata = Glass$Al)
glass_boxcox$Si <- glass_boxcox$Si |> BoxCoxTrans() |> predict(newdata = Glass$Si)
glass_boxcox$Ca <- glass_boxcox$Ca |> BoxCoxTrans() |> predict(newdata = Glass$Ca)glass_boxcox |>
select(-Type) |>
pivot_longer(cols = everything()) |>
ggplot(aes(x = value)) +
geom_histogram(bins = 30, fill = 'skyblue', color = 'black') +
facet_wrap(~ name, scales = 'free')
When viewing the Box-Cox transformed elements, there are notable differences in their center, skew and spread.
Another way to approach the transformation is coded below. The preProcess applies Box-Cox, center, scale, and pca in a pipeline.
glass_trans <- Glass |>
select(-Type) |>
preProcess(method = c('BoxCox', 'center', 'scale', 'pca'))
glass_transCreated from 214 samples and 9 variables
Pre-processing:
- Box-Cox transformation (5)
- centered (9)
- ignored (0)
- principal component signal extraction (9)
- scaled (9)
Lambda estimates for Box-Cox transformation:
-2, -0.1, 0.5, 2, -1.1
PCA needed 7 components to capture 95 percent of the varianceAs before only 5 elements are available for Box-Cox transformation. All were centered, none were ignored, PCA extraction was applied to all, and all were scaled. Additionally, 95% of variance is explained with the first 7 PCA components.
Their loadings are shown below.
glass_trans$rotation PC1 PC2 PC3 PC4 PC5 PC6
RI -0.5273271 0.2931617 -0.18909963 0.13073247 -0.100162337 -0.14493191
Na 0.2330624 0.2943532 0.33954136 0.52437367 0.191509148 0.55703334
Mg -0.1228409 -0.5973021 0.01396475 0.36776189 0.164949703 -0.31381971
Al 0.4390604 0.2389744 -0.30271246 -0.21161203 -0.006701207 0.02714219
Si 0.2027214 -0.1236303 0.53982290 -0.59323075 0.021239416 -0.10240769
K 0.2515777 -0.2330735 -0.59215753 -0.06656522 -0.338650831 0.28413379
Ca -0.4995472 0.3109264 0.04099946 -0.28158708 -0.187405419 0.19102735
Ba 0.2719534 0.4919378 -0.14221977 0.11351848 0.236899299 -0.61638587
Fe -0.1784677 -0.0724442 -0.30521538 -0.28174711 0.848328282 0.24869007
PC7
RI -0.09504494
Na -0.08592573
Mg 0.26083618
Al 0.71816206
Si -0.23388817
K -0.45940175
Ca 0.16140207
Ba -0.31589600
Fe -0.09053409glass_transformed <- predict(glass_trans, Glass |> select(-Type))When plotted as histograms the PCAs are centered around 0, several still show significant skew and outlier values.
glass_transformed |>
pivot_longer(cols = everything()) |>
ggplot(aes(x = value)) +
geom_histogram(bins = 30, fill = 'skyblue', color = 'black') +
facet_wrap(~ name, scales = 'free')
head(glass_transformed) PC1 PC2 PC3 PC4 PC5 PC6 PC7
1 -1.2126444 -0.3942139 0.1730756 1.7193852 -0.1913387 -0.3686905 0.47956412
2 0.6179073 -0.7020476 0.5507034 0.8575350 -0.1566312 0.0618975 0.08525582
3 0.9907027 -0.8876886 0.6452946 0.3027716 -0.1363025 -0.1739302 0.39412058
4 0.1510212 -0.9042336 0.1622361 0.4521567 -0.4291846 -0.2951884 0.10178442
5 0.3582849 -1.0160965 0.5763959 0.1667831 -0.3634192 -0.3289072 -0.13934291
6 0.3408017 -1.3565637 -0.7451275 -1.0568333 1.7762845 0.1433598 0.23899752Correlation
glass_corr <- Glass |>
select(-Type) |>
cor()
glass_corr |> dim()[1] 9 9glass_corr |>
as.data.frame() |>
round(digits = 3) 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.000corrplot(glass_corr, order = 'hclust')
Inspecting the correlations there are several highly correlated combinations of features. The most highly correlated combinations are as follows:
Calcium and Refractive Index are very negatively correlated
Silicon and Refractive Index are very positively correlated
Aluminum and Magnesium are very positively correlated
Barium and Magnesium are very positively correlated.
When a correlation threshold of 0.80 is applied, the correlation between Refractive Index and Calcium is identified as too high; Calcium is recommended for removal from the data based on the chosen threshold.
high_corr <- findCorrelation(glass_corr, cutoff = 0.8)
print(high_corr)[1] 7filtered_glass <- Glass[, -high_corr]
filtered_glass |> as.data.frame() |> head() RI Na Mg Al Si K Ba Fe Type
1 1.52101 13.64 4.49 1.10 71.78 0.06 0 0.00 1
2 1.51761 13.89 3.60 1.36 72.73 0.48 0 0.00 1
3 1.51618 13.53 3.55 1.54 72.99 0.39 0 0.00 1
4 1.51766 13.21 3.69 1.29 72.61 0.57 0 0.00 1
5 1.51742 13.27 3.62 1.24 73.08 0.55 0 0.00 1
6 1.51596 12.79 3.61 1.62 72.97 0.64 0 0.26 1Near Zero Variance
Glass |>
select(-Type) |>
nearZeroVar()integer(0)No Glass features show an issue for zero variance based on the default criteria set in the caret package.
b.
Boxplots & Outlier Analysis
glass_long <- pivot_longer(Glass,
cols = where(is.numeric),
names_to = 'column_name',
values_to = 'value')
ggplot(glass_long, aes(x = column_name, y = value, fill = column_name)) +
geom_boxplot() +
facet_wrap(~ column_name, scales = 'free') +
theme_minimal() +
theme(axis.text.x = element_blank())
Boxplots show that all have observations outside the whiskers, except Mg. When the top and bottom values are examined below though, they indicate continuous percentage values in a reasonable range. When considered inside their context as a percentage makeup of glass, the percentages of the elements all seem plausible.
top_10 <- Glass |>
select(-Type) |>
pivot_longer(
cols = everything(),
names_to = "column",
values_to = "value"
) |>
group_by(column) |>
slice_max(
order_by = value,
n = 10,
with_ties = FALSE
) |>
arrange(desc(value), .by_group = TRUE) |>
mutate(rank = row_number()) |>
ungroup() |>
pivot_wider(
id_cols = rank,
names_from = column,
values_from = value
)
top_10# A tibble: 10 × 10
rank Al Ba Ca Fe K Mg Na RI Si
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 3.5 3.15 16.2 0.51 6.21 4.49 17.4 1.53 75.4
2 2 3.04 2.88 15.0 0.37 6.21 3.98 15.8 1.53 75.2
3 3 3.02 2.2 14.7 0.35 2.7 3.97 15.2 1.53 74.6
4 4 2.88 1.71 14.4 0.34 1.76 3.93 15.0 1.53 74.4
5 5 2.79 1.68 13.4 0.32 1.68 3.9 15.0 1.53 73.9
6 6 2.74 1.67 13.3 0.31 1.46 3.9 15.0 1.53 73.8
7 7 2.68 1.64 13.2 0.3 1.41 3.9 15.0 1.53 73.8
8 8 2.66 1.63 12.5 0.29 1.1 3.89 14.9 1.53 73.7
9 9 2.54 1.59 12.2 0.28 0.97 3.87 14.9 1.52 73.7
10 10 2.51 1.59 11.6 0.28 0.81 3.86 14.9 1.52 73.6bot_10 <- Glass |>
select(-Type) |>
pivot_longer(
cols = everything(),
names_to = "column",
values_to = "value"
) |>
group_by(column) |>
slice_min(
order_by = value,
n = 10,
with_ties = FALSE
) |>
arrange(value, .by_group = TRUE) |>
mutate(rank = row_number()) |>
ungroup() |>
pivot_wider(
id_cols = rank,
names_from = column,
values_from = value
)
bot_10# A tibble: 10 × 10
rank Al Ba Ca Fe K Mg Na RI Si
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.29 0 5.43 0 0 0 10.7 1.51 69.8
2 2 0.34 0 5.79 0 0 0 11.0 1.51 69.9
3 3 0.47 0 5.87 0 0 0 11.0 1.51 70.2
4 4 0.47 0 6.47 0 0 0 11.2 1.51 70.3
5 5 0.51 0 6.65 0 0 0 11.4 1.51 70.4
6 6 0.56 0 6.93 0 0 0 11.6 1.51 70.5
7 7 0.56 0 6.96 0 0 0 12.0 1.51 70.6
8 8 0.58 0 7.08 0 0 0 12.2 1.52 70.7
9 9 0.65 0 7.36 0 0 0 12.2 1.52 71.2
10 10 0.66 0 7.59 0 0 0 12.3 1.52 71.2c.
As indicated above, there are several options for transformation. Considering the quantity of zero values, a log based transformation is not practical with this data. Additionally, the zero values do not indicate missing data, but are indicators of that element not existing as part of the glass.
That being said, Box-Cox is available for 5 of the 9 elements, and centering and scaling and PCA all appear viable options for transformation, resulting in 7 PCs capable of explaining 95% of the variation.
Exercise 3.2
a.
data(Soybean)str(Soybean)'data.frame': 683 obs. of 36 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 ...
$ crop.hist : Factor w/ 4 levels "0","1","2","3": 2 3 2 2 3 4 3 2 4 3 ...
$ area.dam : Factor w/ 4 levels "0","1","2","3": 2 1 1 1 1 1 1 1 1 1 ...
$ sever : Factor w/ 3 levels "0","1","2": 2 3 3 3 2 2 2 2 2 3 ...
$ seed.tmt : Factor w/ 3 levels "0","1","2": 1 2 2 1 1 1 2 1 2 1 ...
$ germ : Ord.factor w/ 3 levels "0"<"1"<"2": 1 2 3 2 3 2 1 3 2 3 ...
$ plant.growth : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
$ leaves : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
$ leaf.halo : Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 1 1 1 1 ...
$ leaf.marg : Factor w/ 3 levels "0","1","2": 3 3 3 3 3 3 3 3 3 3 ...
$ leaf.size : Ord.factor w/ 3 levels "0"<"1"<"2": 3 3 3 3 3 3 3 3 3 3 ...
$ leaf.shread : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ leaf.malf : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ leaf.mild : Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 1 1 1 1 ...
$ stem : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
$ lodging : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 2 1 1 1 ...
$ stem.cankers : Factor w/ 4 levels "0","1","2","3": 4 4 4 4 4 4 4 4 4 4 ...
$ canker.lesion : Factor w/ 4 levels "0","1","2","3": 2 2 1 1 2 1 2 2 2 2 ...
$ fruiting.bodies: Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
$ ext.decay : Factor w/ 3 levels "0","1","2": 2 2 2 2 2 2 2 2 2 2 ...
$ mycelium : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ int.discolor : Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 1 1 1 1 ...
$ sclerotia : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ fruit.pods : Factor w/ 4 levels "0","1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
$ fruit.spots : Factor w/ 4 levels "0","1","2","4": 4 4 4 4 4 4 4 4 4 4 ...
$ seed : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ mold.growth : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ seed.discolor : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ seed.size : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ shriveling : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ roots : Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 1 1 1 1 ...library(janitor)
frequency_table <- Soybean |>
select(-Class) |>
lapply(tabyl) |> adorn_pct_formatting()
head(frequency_table)$date
X[[i]] n percent valid_percent
0 26 3.8% 3.8%
1 75 11.0% 11.0%
2 93 13.6% 13.6%
3 118 17.3% 17.3%
4 131 19.2% 19.2%
5 149 21.8% 21.8%
6 90 13.2% 13.2%
<NA> 1 0.1% -
$plant.stand
X[[i]] n percent valid_percent
0 354 51.8% 54.7%
1 293 42.9% 45.3%
<NA> 36 5.3% -
$precip
X[[i]] n percent valid_percent
0 74 10.8% 11.5%
1 112 16.4% 17.4%
2 459 67.2% 71.2%
<NA> 38 5.6% -
$temp
X[[i]] n percent valid_percent
0 80 11.7% 12.3%
1 374 54.8% 57.3%
2 199 29.1% 30.5%
<NA> 30 4.4% -
$hail
X[[i]] n percent valid_percent
0 435 63.7% 77.4%
1 127 18.6% 22.6%
<NA> 121 17.7% -
$crop.hist
X[[i]] n percent valid_percent
0 65 9.5% 9.7%
1 165 24.2% 24.7%
2 219 32.1% 32.8%
3 218 31.9% 32.7%
<NA> 16 2.3% -Inspecting the frequency distributions produced above, there are many features which show a possible near-zero variance. Some low variation features include leaf.mild, mycelium, sclerotia, seed size, shriveling, and more.
Soybean |>
select(-Class) |>
nearZeroVar()[1] 18 25 27soybean_no_class <- Soybean |>
select(-Class)
head(soybean_no_class[,c(18,25,27)]) leaf.mild mycelium sclerotia
1 0 0 0
2 0 0 0
3 0 0 0
4 0 0 0
5 0 0 0
6 0 0 0As suspected upon the review of the data, the nearZeroVar function identified leaf.mild, mycelium, and sclerotia as all being near-zero variance.
b.
Missing Values
total_pieces_of_data <- prod(dim(soybean_no_class))
total_na_values <- sum(is.na(soybean_no_class))
print(total_pieces_of_data)[1] 23905print(total_na_values)[1] 2337print(total_na_values / total_pieces_of_data)[1] 0.09776197About 9.8% of all feature data is marked as NA.
soybean_no_class |>
summarise_all(~ mean(is.na(.)) * 100) |>
pivot_longer(cols = everything(), names_to = 'feature', values_to = 'pct_na') |>
arrange(desc(pct_na))# A tibble: 35 × 2
feature pct_na
<chr> <dbl>
1 hail 17.7
2 sever 17.7
3 seed.tmt 17.7
4 lodging 17.7
5 germ 16.4
6 leaf.mild 15.8
7 fruiting.bodies 15.5
8 fruit.spots 15.5
9 seed.discolor 15.5
10 shriveling 15.5
# ℹ 25 more rowsThe percentage of missing data for each feature ranges from near 18% to 0% missing. There is not one or even a few columns that have execess amounts of missing data. Although many have near identical rates of missing data. This might imply a correlary reason that those features to not be recorded. Further information about how the data is collected would help clarify why some have identical missing rates.
na_sum <- Soybean |>
group_by(Class) |>
summarise(across(everything(), ~ sum(is.na(.)))) |>
mutate(total_na = rowSums(pick(-Class))) |>
select(Class, total_na) |>
filter(total_na != 0)
na_sum# A tibble: 5 × 2
Class total_na
<fct> <dbl>
1 2-4-d-injury 450
2 cyst-nematode 336
3 diaporthe-pod-&-stem-blight 177
4 herbicide-injury 160
5 phytophthora-rot 1214n_distinct(Soybean$Class)[1] 19An examination of the NA values indicates that they only exist for five of the nineteen classes. All NA values are only in five classes.
c.
Inspecting the data set, the classes appear to be disease or damage of soybean plants. It would be valuable to know more about each type of disease or injury. Considering that the NA values are relegated to just five classes, there is the potential that the NA data may not be possible to collect for those specific classes.
Furthermore, since the NA values often take up the majority of entries for that class type within a feature, imputation may be misleading. NA values are also spread across all features except one, so eliminating features is not practical. It would be possible to set a maximum threshold of NA values and eliminate any feature that exceeds that level, but ultimately this might end up in the excisement of useful predictors.
In order to deal with this missing data there is a third procedure that is possible:
Code NA values as their own value which indicates the missing nature of the value, coding “missingness” as its own category.
If the feature uses ordered numerical indicators, impute and appropriate value, paired with a missing category.
Considering the NA values appear for only 5 Classes, the fact that NA was recorded may be itself a useful predictor, so finding a way to retain this information could be valuable.