pkgs <- c("mlbench", "caret", "e1071")
new <- pkgs[!(pkgs %in% installed.packages()[, "Package"])]
if (length(new)) install.packages(new)Data624Homework4
3.1
# Loading required packages
library(mlbench)
library(caret)
library(e1071)
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
par(mfrow = c(3, 3), mar = c(3, 3, 2, 1))
for (v in names(Glass)[1:9]) {
hist(Glass[[v]], main = v, xlab = "", col = "blue")
}
par(mfrow = c(1, 1))The histograms of the 9 predictors show that NA and Si are roughly symetric, while Ri, Ai, and Ca are kind of skewed. Mg is bimodal the other groups of glasses are at 0. The spread of the x axes for K, Ba, and Ca might point at outliers.
pairs(Glass[, 1:9], pch = 20, cex = 0.5)
round(cor(Glass[, 1:9]), 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.00The scatterplot shows that RI and Ca have strong positive linear relationship. Ri has a weaker negative relationship with Si. The other parts shoows no clear pattern of collinearity overall.
b
par(mfrow = c(3, 3), mar = c(3, 3, 2, 1))
for (v in names(Glass)[1:9]) {
boxplot(Glass[[v]], main = v, horizontal = TRUE)
}
par(mfrow = c(1, 1))This boxplot show potential outliers in all predictors except Mg
c
sk <- apply(Glass[, 1:9], 2, skewness)
round(sk, 2) RI Na Mg Al Si K Ca Ba Fe
1.60 0.45 -1.14 0.89 -0.72 6.46 2.02 3.37 1.73 These results agree with the histograms and the outliers seen in the boxplots
d
colSums(Glass[, 1:9] == 0) RI Na Mg Al Si K Ca Ba Fe
0 0 42 0 0 30 0 176 144 Because of the zeros on 4, a box box transformation cannot be applied to them. but for those that are positive, those are candidates for Box Cox
pos <- c("RI", "Na", "Al", "Si", "Ca")
lambdas <- sapply(pos, function(v) BoxCoxTrans(Glass[[v]])$lambda)
round(lambdas, 2) RI Na Al Si Ca
-2.0 -0.1 0.5 2.0 -1.1 bc <- sapply(pos, function(v) {
predict(BoxCoxTrans(Glass[[v]]), Glass[[v]])
})
round(rbind(before = sk[pos], after = apply(bc, 2, skewness)), 2) RI Na Al Si Ca
before 1.60 0.45 0.89 -0.72 2.02
after 1.57 0.03 0.09 -0.65 -0.19Applying Box Cox reduced skewness substantially for Na (0.45 to 0.03), Al (0.89 to 0.09), and Ca (2.02 to negative 0.19). RI (1.60 to 1.57) and Si (negative 0.72 to negative 0.65) were essentially unchanged, because their estimated lambdas hit the limits of the search grid and their skew is driven by a few extreme observations. For these, a spatial sign transformation, which handles outliers, may be more appropriate.
ss <- spatialSign(scale(Glass[, 1:9]))
head(round(ss, 3), 3) RI Na Mg Al Si K Ca Ba Fe
1 0.386 0.126 0.555 -0.306 -0.499 -0.297 -0.064 -0.156 -0.259
2 -0.178 0.423 0.455 -0.122 0.073 -0.019 -0.568 -0.252 -0.419
3 -0.474 0.099 0.395 0.125 0.288 -0.108 -0.545 -0.232 -0.385A spatial sign transformation was applied after centering and scaling all the 9 predictors, which projects every observation onto a unit sphere. This limits the influence of extreme points, since an observation with large values is divided by its own large length. Unlike Box Cox, it can be applied to Mg, K, Ba, and Fe despite their zeros, and it is a reasonable choice for RI and Si, where Box Cox barely reduced skewness.
range(sqrt(rowSums(ss^2)))[1] 1 13.2
data(Soybean)
dim(Soybean)[1] 683 36str(Soybean[, 1:6])'data.frame': 683 obs. of 6 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 ...a
tabs <- lapply(Soybean[, -1], table)
tabs[1:3]$date
0 1 2 3 4 5 6
26 75 93 118 131 149 90
$plant.stand
0 1
354 293
$precip
0 1 2
74 112 459 Frequency tables for the first three predictors show no degenerate variables. date is spread across seven levels, plant.stand is nearly balanced and precip is skewed.
top <- sapply(Soybean[, -1], function(x) max(table(x)) / sum(!is.na(x)))
round(sort(top, decreasing = TRUE), 2) mycelium sclerotia shriveling leaf.mild lodging
0.99 0.97 0.93 0.93 0.93
leaf.malf int.discolor seed.size seed.discolor leaves
0.92 0.90 0.90 0.89 0.89
mold.growth roots leaf.shread fruiting.bodies seed
0.89 0.85 0.84 0.82 0.81
hail ext.decay precip fruit.pods plant.growth
0.77 0.77 0.71 0.68 0.66
fruit.spots leaf.marg stem.cankers sever temp
0.60 0.60 0.59 0.57 0.57
leaf.halo stem plant.stand leaf.size seed.tmt
0.57 0.56 0.55 0.55 0.54
canker.lesion germ area.dam crop.hist date
0.50 0.37 0.33 0.33 0.22 nzv <- nearZeroVar(Soybean[, -1], saveMetrics = TRUE)
nzv[nzv$nzv, ] freqRatio percentUnique zeroVar nzv
leaf.mild 26.75 0.4392387 FALSE TRUE
mycelium 106.50 0.2928258 FALSE TRUE
sclerotia 31.25 0.2928258 FALSE TRUEUsing carets nearZeroVar with the standard cutoffs the three Soybean predictors are degenerate. Each is dominated by a single level, so it provides almost no information for separating the 19 classes and is a candidate for removal before modeling.
na_count <- colSums(is.na(Soybean[, -1]))
sort(na_count, decreasing = TRUE) hail sever seed.tmt lodging germ
121 121 121 121 112
leaf.mild fruiting.bodies fruit.spots seed.discolor shriveling
108 106 106 106 106
leaf.shread seed mold.growth seed.size leaf.halo
100 92 92 92 84
leaf.marg leaf.size leaf.malf fruit.pods precip
84 84 84 84 38
stem.cankers canker.lesion ext.decay mycelium int.discolor
38 38 38 38 38
sclerotia plant.stand roots temp crop.hist
38 36 31 30 16
plant.growth stem date area.dam leaves
16 16 1 1 0 Missing values are widespread in the Soybean predictors, many predictors share exactly the same number of missing values which suggests that groups of predictors tend to be missing together on the same observations rather than at random. Only date, area.dam, and leaves are nearly complete.
any_na <- rowSums(is.na(Soybean[, -1])) > 0
miss_by_class <- tapply(any_na, Soybean$Class, mean)
round(sort(miss_by_class, decreasing = TRUE), 2) 2-4-d-injury cyst-nematode
1.00 1.00
diaporthe-pod-&-stem-blight herbicide-injury
1.00 1.00
phytophthora-rot alternarialeaf-spot
0.77 0.00
anthracnose bacterial-blight
0.00 0.00
bacterial-pustule brown-spot
0.00 0.00
brown-stem-rot charcoal-rot
0.00 0.00
diaporthe-stem-canker downy-mildew
0.00 0.00
frog-eye-leaf-spot phyllosticta-leaf-spot
0.00 0.00
powdery-mildew purple-seed-stain
0.00 0.00
rhizoctonia-root-rot
0.00 table(Soybean$Class[any_na])
2-4-d-injury alternarialeaf-spot
16 0
anthracnose bacterial-blight
0 0
bacterial-pustule brown-spot
0 0
brown-stem-rot charcoal-rot
0 0
cyst-nematode diaporthe-pod-&-stem-blight
14 15
diaporthe-stem-canker downy-mildew
0 0
frog-eye-leaf-spot herbicide-injury
0 8
phyllosticta-leaf-spot phytophthora-rot
0 68
powdery-mildew purple-seed-stain
0 0
rhizoctonia-root-rot
0 sum(any_na)[1] 121The first four classes are missing values in every observation, so deleting incomplete rows would remove them from the data entirely and leave phytophthora-rot with only about 20 observations.
cc <- Soybean[complete.cases(Soybean), ]
cbind(all = table(Soybean$Class), complete = table(cc$Class)) all complete
2-4-d-injury 16 0
alternarialeaf-spot 91 91
anthracnose 44 44
bacterial-blight 20 20
bacterial-pustule 20 20
brown-spot 92 92
brown-stem-rot 44 44
charcoal-rot 20 20
cyst-nematode 14 0
diaporthe-pod-&-stem-blight 15 0
diaporthe-stem-canker 20 20
downy-mildew 20 20
frog-eye-leaf-spot 91 91
herbicide-injury 8 0
phyllosticta-leaf-spot 20 20
phytophthora-rot 88 20
powdery-mildew 20 20
purple-seed-stain 20 20
rhizoctonia-root-rot 20 20A model trained on the complete cases could never predict these diseases, row deletion is not an acceptable strategy for this dataset.
drop_cols <- c("mycelium", "sclerotia", "leaf.mild")
soy <- Soybean[, !(names(Soybean) %in% drop_cols)]
soy[, -1] <- lapply(soy[, -1], function(x) {
x <- addNA(x)
levels(x)[is.na(levels(x))] <- "missing"
x
})
colSums(is.na(soy))[1:5] Class date plant.stand precip temp
0 0 0 0 0 After recoding, no missing values remain in any predictor, so all 683 observations, including the four classes that would otherwise be lost, are kept for modeling
table(soy$hail)
0 1 missing
435 127 121 In the Soybean data, caret’s nearZeroVar flagged 3 degenerate predictors: mycelium, sclerotia, and leaf.mild. Missing values are strongly tied to the class, since 4 classes have them in every observation and 14 classes have none.
Deleting incomplete rows would remove four classes entirely, so I dropped the 3 degenerate predictors and recoded missing values as their own factor level, which keeps all 683 observations.