Классификация
Небожатко Екатерина
library(MASS)
library(psych)
library(klaR)
library(ICSNP)
library(mvnormtest)
library(HDMD)
library(biotools)## ---
## biotools version 3.0library(Hotelling)
library(candisc)Данные
df <- read.csv("/home/ekaterina/Downloads/winequality-red_good.csv")
#df$quality <- as.factor(df$quality)
df1 <- df[,-c(1,3,5,6,8)]
head(df1)## volatile.acidity residual.sugar total.sulfur.dioxide pH sulphates
## 1 0.70 1.9 34 3.51 0.56
## 2 0.88 2.6 67 3.20 0.68
## 3 0.76 2.3 54 3.26 0.65
## 4 0.28 1.9 60 3.16 0.58
## 5 0.70 1.9 34 3.51 0.56
## 6 0.66 1.8 40 3.51 0.56
## alcohol quality
## 1 9.4 5
## 2 9.8 5
## 3 9.8 5
## 4 9.8 6
## 5 9.4 5
## 6 9.4 5Описание признаков
residual sugar – остаточный сахар. Чем меньше, тем лучше, это означает, что скорее всего не добавляли сахарных сиропов и прочего, ведь такого в вине не должно быть в принципе
pH – показатель кислотности
total sulfur dioxide – полное содержание оксида серы (поступает из винограда и добавляют в качестве пищевой добавки-консерванта)
sulphates – содержание сульфатов (минеральная соль) – чем больше, тем лучше
volatile acidity – летучая кислотность (кислотность летучих кислот) – чем меньше, тем лучше, ибо повышенное содержание говорит о том, что вино скисало или имеет место быть наличие бактерий
quality - качество вина по шкале от 3 до 8
pairs.panels(df1, lm=TRUE, bg=c("red","yellow","blue", "pink", "green")[df1$quality], pch=21, lwd = 0)
Вин среднего качества больше всего и они перекрывают вины более низкого и высокого качества.
Прологорифмируем некоторые переменные
df1$volatile.acidity <- log(df1$volatile.acidity)
df1$residual.sugar <- log(df1$residual.sugar)
df1$sulphates <- log(df1$sulphates)
df1$total.sulfur.dioxide <- log(df1$total.sulfur.dioxide)
df1$alcohol <- log(df1$alcohol)
pairs.panels(df1, lm=TRUE, bg=c("red","yellow","blue", "pink", "green")[df1$quality], pch=21, lwd = 0)
Проверка модели
Проверка на нормальность
normTest <- function (x){
res <- shapiro.test(x)
return(c(res$statistic, p.value = res$p.value)) }
t(sapply(df1, normTest))## W p.value
## volatile.acidity 0.9864417 4.210479e-11
## residual.sugar 0.8550736 5.471063e-36
## total.sulfur.dioxide 0.9899255 4.632268e-09
## pH 0.9934863 1.712237e-06
## sulphates 0.9588877 7.623993e-21
## alcohol 0.9464349 1.155071e-23
## quality 0.8575895 9.515085e-36Гипотезы о нормальном распределении отвергаются для всех признаков. Однако, это не так страшно, потому что все выводы остаются верны асимптотически.
Гомоскедастичность
quality <- df1[,7]
Y <- df1[,-7]
boxM(Y, quality)##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: Y
## Chi-Sq (approx.) = 294.48, df = 105, p-value < 2.2e-16Гипотеза о гомоскедастичности отвергается.
Расстояние Махаланобиса между группами
quality1 <- as.factor(quality)
md <- pairwise.mahalanobis(Y, quality1)$distance
colnames(md) <- levels(quality1)
rownames(md) <- levels(quality1)
print(md, digits = 3)## 3 4 5 6 7 8
## 3 0.000 0.171 0.8939 0.6755 0.7597 0.6526
## 4 0.171 0.000 0.3006 0.1821 0.3161 0.2609
## 5 0.894 0.301 0.0000 0.0995 0.4035 0.3814
## 6 0.675 0.182 0.0995 0.0000 0.1099 0.0989
## 7 0.760 0.316 0.4035 0.1099 0.0000 0.0196
## 8 0.653 0.261 0.3814 0.0989 0.0196 0.0000Расстояние между 5 и 6, 6 и 7, 6 и 8, 7 и 8 очень маленькое, почти ноль. То есть вины выше среднего почти не различаются по расстояние Махалонобиса, чего не наблюдается у вин ниже среднего. Инетересно, что наибольшее расстояние между 3 и 5 качеством, а не между 3 и 8, как можно было этого ожидать.
MANOVA
fit <- manova(as.matrix(Y)~quality)
summary(fit, test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## quality 1 0.65108 142.2 6 1592 < 2.2e-16 ***
## Residuals 1597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Гипотеза о значимости дискриминантного анализа не отвергается
Дискриминантный анализ
LDA
Предполагаем гомоскедастичность, хотя гипотеза у нас отверглась
Везде будем использовать априорные вероятности пропорциональные количеству индивидов в группе, потому что у нас группы очень сильно различаются по объему. Они идут по умолчанию в методе lda
#### LDA ####
wine.lda <- lda(Y, quality)
wine.lda## Call:
## lda(Y, quality)
##
## Prior probabilities of groups:
## 3 4 5 6 7 8
## 0.006253909 0.033145716 0.425891182 0.398999375 0.124452783 0.011257036
##
## Group means:
## volatile.acidity residual.sugar total.sulfur.dioxide pH sulphates
## 3 -0.1845772 0.8601003 3.011119 3.398000 -0.5807948
## 4 -0.4197593 0.8743572 3.309623 3.381509 -0.5625070
## 5 -0.5920088 0.8489759 3.795272 3.304949 -0.5049719
## 6 -0.7535112 0.8295904 3.528611 3.318072 -0.4148078
## 7 -0.9680623 0.9127478 3.302340 3.290754 -0.3160045
## 8 -0.9051405 0.8616992 3.279611 3.267222 -0.2739620
## alcohol
## 3 2.294936
## 4 2.324793
## 5 2.289943
## 6 2.358881
## 7 2.435847
## 8 2.487785
##
## Coefficients of linear discriminants:
## LD1 LD2 LD3 LD4
## volatile.acidity -1.18305241 -1.73817511 1.3454262 1.8332402
## residual.sugar 0.09120509 -0.38436959 1.3136576 -2.1942244
## total.sulfur.dioxide -0.33302480 1.04499204 0.1629951 0.5860921
## pH -0.62480306 -1.95217133 -5.0961882 -2.0084859
## sulphates 1.72875250 -0.02314187 -0.8756144 1.0528026
## alcohol 8.80996913 -1.41198452 4.3805380 5.4975821
## LD5
## volatile.acidity 0.60167850
## residual.sugar -0.01722428
## total.sulfur.dioxide -0.67846682
## pH -1.99964105
## sulphates 3.39886384
## alcohol -5.30074204
##
## Proportion of trace:
## LD1 LD2 LD3 LD4 LD5
## 0.9094 0.0726 0.0101 0.0066 0.0014Качество классификации
wine.ldap <- predict(wine.lda, Y)$class
t1 <- table(wine.ldap, quality) #classification quality
t1## quality
## wine.ldap 3 4 5 6 7 8
## 3 1 0 0 0 0 0
## 4 2 1 2 0 0 0
## 5 6 34 512 229 11 0
## 6 1 17 161 345 107 9
## 7 0 1 6 64 81 9
## 8 0 0 0 0 0 0Доля ошибок
1-signif(diag(prop.table(t1, 2)),3)## 3 4 5 6 7 8
## 0.9000 0.9811 0.2480 0.4590 0.5930 1.0000Качество классификации оставляет желать лучшего. В 8 группу не поплал никто, однако полностью правильно определена 3 группа, хоть туда и попали индивиды из других групп. Больше всего ошибок у 4 группы. Она оказалась распределена между 3 и 5.
LDA Cross-validation
wine.lda.cv <- lda(Y, quality, CV = TRUE)
t2 <- table(wine.lda.cv$class, quality)
t2## quality
## 3 4 5 6 7 8
## 3 0 0 0 0 0 0
## 4 3 0 2 0 0 0
## 5 6 34 510 230 11 0
## 6 1 18 163 342 110 9
## 7 0 1 6 66 78 9
## 8 0 0 0 0 0 0Доля ошибок
1-signif(diag(prop.table(t2, 2)), 3)## 3 4 5 6 7 8
## 1.000 1.000 0.251 0.464 0.608 1.000Проверка по cross-validation показала совсем плохие результаты. 3 и 8 класс отсутствуют, в 4 никто правильно не попал.
QDA
Используем теперь предположение о гетероскедастичности модели.
##### QDA #####
wine.qda <- qda(Y, quality)
wine.qda## Call:
## qda(Y, quality)
##
## Prior probabilities of groups:
## 3 4 5 6 7 8
## 0.006253909 0.033145716 0.425891182 0.398999375 0.124452783 0.011257036
##
## Group means:
## volatile.acidity residual.sugar total.sulfur.dioxide pH sulphates
## 3 -0.1845772 0.8601003 3.011119 3.398000 -0.5807948
## 4 -0.4197593 0.8743572 3.309623 3.381509 -0.5625070
## 5 -0.5920088 0.8489759 3.795272 3.304949 -0.5049719
## 6 -0.7535112 0.8295904 3.528611 3.318072 -0.4148078
## 7 -0.9680623 0.9127478 3.302340 3.290754 -0.3160045
## 8 -0.9051405 0.8616992 3.279611 3.267222 -0.2739620
## alcohol
## 3 2.294936
## 4 2.324793
## 5 2.289943
## 6 2.358881
## 7 2.435847
## 8 2.487785Качество классификации
wine.qdap <- predict(wine.qda, Y)$class
t3 <- table(wine.qdap, quality) #classification quality
t3## quality
## wine.qdap 3 4 5 6 7 8
## 3 4 0 4 1 0 0
## 4 0 5 5 7 0 0
## 5 5 33 533 246 18 0
## 6 1 14 130 343 102 11
## 7 0 1 9 41 75 4
## 8 0 0 0 0 4 3Доля ошибок
1 - signif(diag(prop.table(t3, 2)),3)## 3 4 5 6 7 8
## 0.6000 0.9057 0.2170 0.4620 0.6230 0.8330Качество классификаии радует чуть больше, потому что у нас хотя бы появились 3 и 8 классы. Однако, количество верно определенных индивидов все еще очень небольшое.
QDA Cross-validation
wine.qda.cv <- qda(Y, quality, CV = TRUE)
t4 <- table(wine.qda.cv$class, quality)
t4## quality
## 3 4 5 6 7 8
## 3 0 0 5 1 0 0
## 4 2 1 5 7 0 0
## 5 7 36 527 252 18 0
## 6 1 14 135 334 106 11
## 7 0 2 9 44 69 6
## 8 0 0 0 0 6 1Доля ошибок
1 - signif(diag(prop.table(t4, 2)), 3)## 3 4 5 6 7 8
## 1.0000 0.9811 0.2260 0.4760 0.6530 0.9444quality.factor <- as.factor(quality)Проверка по cross-validation практически стерла 3 и 4 классы.
Делаем вывод, что проблемные классы - это наиболее плохое и хорошее качество вина (3,4 и 8). Принадлежность к остальным определяет примерно одниково для всех методов и проверок и составляет около 50% верных попаданий.
LDA, out-of-sample
Разделим теперь выборку на обучающую и тестируемую.
#### LDA, out-of-sample ####
wine.train <- df1[-seq(1, nrow(df1), 5),]
wine.unknown <- df1[seq(1,nrow(df1),5),]
wine.lda2 <- lda(wine.train[,1:6], wine.train[,7])
wine.lda2p <- predict(wine.lda2, wine.unknown[,1:6])$class
t5 <- table(wine.lda2p, wine.unknown[,7])
signif(prop.table(t5,1), 3)##
## wine.lda2p 3 4 5 6 7 8
## 3
## 4 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## 5 0.0061 0.0183 0.6950 0.2800 0.0000 0.0000
## 6 0.0000 0.0161 0.2020 0.5560 0.1940 0.0323
## 7 0.0000 0.0323 0.0000 0.4520 0.4520 0.0645
## 8signif(diag(prop.table(t5,1)), 3)## 3 4 5 6 7 8
## NaN 0.000 0.695 0.556 0.452 NaNКачество классификации не улучшилось.
Объединение групп
Объединим группы по расстоянию Махаланобиса. Пусть их будет 3:
1 - низкое качество 3 и 4
2 - среднее 5 и 6
3 - высокое 7 и 8
df1$new_quality <- as.factor(df1$quality)
levels(df1$new_quality) <- c(1,1,2,2,3,3)
#df1$new_quality <- as.numeric(df1$new_quality)
quality2 <- df1[,8]
pairs.panels(df1[,-7], bg=c("red","yellow","blue")[df1$new_quality], pch=21, lwd = 0)
Количество индивидов по группам
sapply(levels(df1$new_quality), function(level){
nrow(df1[df1$new_quality == level,])
})## 1 2 3
## 63 1319 217Канонический анализ
MANOVA
manova_wine1 <- manova(as.matrix(df1[-c(7,8)])~df1$new_quality, data = df1)
summary(manova_wine1, test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## df1$new_quality 2 0.72537 46.177 12 3182 < 2.2e-16 ***
## Residuals 1596
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Дискриминантный анализ значим
Избыточность признаков
greedy.wilks(df1[-c(7,8)], quality2)## Formula containing included variables:
##
## quality2 ~ alcohol + volatile.acidity + sulphates + total.sulfur.dioxide +
## pH + residual.sugar
## <environment: 0x8397120>
##
##
## Values calculated in each step of the selection procedure:
##
## vars Wilks.lambda F.statistics.overall p.value.overall
## 1 alcohol 0.8356592 156.93477 6.015734e-63
## 2 volatile.acidity 0.7702488 111.18878 7.601321e-89
## 3 sulphates 0.7515643 81.55872 3.173240e-95
## 4 total.sulfur.dioxide 0.7353268 66.17501 8.871697e-101
## 5 pH 0.7283892 54.67085 3.482194e-102
## 6 residual.sugar 0.7253675 46.17689 7.430611e-102
## F.statistics.diff p.value.diff
## 1 156.934772 6.015734e-63
## 2 67.724532 0.000000e+00
## 3 19.814099 3.162325e-09
## 4 17.588387 2.783078e-08
## 5 7.581537 5.283871e-04
## 6 3.313825 3.662787e-02Не значим residual.sugar. Уберем его и посмотрим, как изменилась статистика
MANOVA еще раз
manova_wine2 <- manova(as.matrix(df1[-c(2,7,8)])~df1$new_quality, data = df1)
summary(manova_wine2, test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## df1$new_quality 2 0.72839 54.671 10 3184 < 2.2e-16 ***
## Residuals 1596
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Чуть-чуть увеличилась статистика. Осталось значимо
Hotelling test
df2 <- df1[,-c(2,7)]
#Между 1 и 2 группой
print(hotelling.test(.~new_quality, data = df2, pair = c(1,2)))## Test stat: 15.454
## Numerator df: 5
## Denominator df: 1376
## P-value: 8.216e-15#Между 1 и 3 группой
print(hotelling.test(.~new_quality, data = df2, pair = c(1,3)))## Test stat: 62.59
## Numerator df: 5
## Denominator df: 274
## P-value: 0#Между 2 и 3 группой
print(hotelling.test(.~new_quality, data = df2, pair = c(2,3)))## Test stat: 95.027
## Numerator df: 5
## Denominator df: 1530
## P-value: 0Все группы значимо различаются между собой
CDA
cd <- candisc(manova_wine2, data = df1[,-c(2,7)])
summary(cd)##
## Canonical Discriminant Analysis for df1$new_quality:
##
## CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.242718 0.320511 0.28084 88.987 88.987
## 2 0.038154 0.039667 0.28084 11.013 100.000
##
## Class means:
##
## Can1 Can2
## 1 -0.72635 -0.948692
## 2 -0.19623 0.060328
## 3 1.40360 -0.091266
##
## std coefficients:
## Can1 Can2
## volatile.acidity -0.38640 -0.574567
## total.sulfur.dioxide -0.12495 0.704223
## pH -0.19536 -0.204733
## sulphates 0.31366 0.028481
## alcohol 0.73549 -0.257746cd$coeffs.std## Can1 Can2
## volatile.acidity -0.3864049 -0.57456671
## total.sulfur.dioxide -0.1249535 0.70422272
## pH -0.1953590 -0.20473323
## sulphates 0.3136629 0.02848083
## alcohol 0.7354916 -0.25774632cd$structure## Can1 Can2
## volatile.acidity -0.6494908 -0.5187061
## total.sulfur.dioxide -0.2957879 0.7163536
## pH -0.1480944 -0.4073705
## sulphates 0.5030593 0.2354997
## alcohol 0.8113886 -0.3452305print(cd)##
## Canonical Discriminant Analysis for df1$new_quality:
##
## CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.242718 0.320511 0.28084 88.987 88.987
## 2 0.038154 0.039667 0.28084 11.013 100.000
##
## Test of H0: The canonical correlations in the
## current row and all that follow are zero
##
## LR test stat approx F numDF denDF Pr(> F)
## 1 0.72839 54.671 10 3184 < 2.2e-16 ***
## 2 0.96185 15.797 4 1593 1.098e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Все значимо
plot(cd)## Vector scale factor set to 4.586
На графике видно, как сильно перемешана 1 группа (зеленый цвет) со всеми остальными. Это может объяснить очень плохую классификацию по 1 группе.
Классификация. Априорные вероятности пропорциональны размерам групп
Сравним между собой только QDA и LDA. Будем использовать априорные вероятности, пропорциональные объему групп.
LDA
Y2 <- df1[-c(2,7,8)]
wine.lda <- lda(Y2, quality2)
wine.lda## Call:
## lda(Y2, quality2)
##
## Prior probabilities of groups:
## 1 2 3
## 0.03939962 0.82489056 0.13570982
##
## Group means:
## volatile.acidity total.sulfur.dioxide pH sulphates alcohol
## 1 -0.3824288 3.262241 3.384127 -0.5654098 2.320054
## 2 -0.6701275 3.666288 3.311296 -0.4613595 2.323288
## 3 -0.9628430 3.300454 3.288802 -0.3125171 2.440156
##
## Coefficients of linear discriminants:
## LD1 LD2
## volatile.acidity -1.1590714 1.7234871
## total.sulfur.dioxide -0.1807566 -1.0187218
## pH -1.2720321 1.3330697
## sulphates 1.4489051 -0.1315617
## alcohol 8.1301265 2.8491288
##
## Proportion of trace:
## LD1 LD2
## 0.8899 0.1101wine.ldap <- predict(wine.lda, Y2)$class
t1 <- table(wine.ldap, quality2)
t1## quality2
## wine.ldap 1 2 3
## 1 1 2 0
## 2 62 1254 133
## 3 0 63 84Доля ошибок
1-signif(diag(prop.table(t1, 2)),3)## 1 2 3
## 0.9841 0.0490 0.6130На матриксплоте видно насколько смешаны данные. Было бы странно ожидать каких-то иных результатов. Очень многие индивиды были ошибочно приписаны ко 2 группе.
Попарный дискриминантный анализ.
partimat(new_quality ~ ., df2, method="lda")
Здесь видно, что по некоторым признакам 1 группа вообще не представлена.
QDA
wine.qda <- qda(Y2, quality2)
wine.qda## Call:
## qda(Y2, quality2)
##
## Prior probabilities of groups:
## 1 2 3
## 0.03939962 0.82489056 0.13570982
##
## Group means:
## volatile.acidity total.sulfur.dioxide pH sulphates alcohol
## 1 -0.3824288 3.262241 3.384127 -0.5654098 2.320054
## 2 -0.6701275 3.666288 3.311296 -0.4613595 2.323288
## 3 -0.9628430 3.300454 3.288802 -0.3125171 2.440156wine.qdap <- predict(wine.qda, Y2)$class
t2 <- table(wine.qdap, quality2)
t2## quality2
## wine.qdap 1 2 3
## 1 6 9 0
## 2 57 1237 125
## 3 0 73 92Доля ошибок
1-signif(diag(prop.table(t2, 2)),3)## 1 2 3
## 0.9048 0.0620 0.5760Доля ошибок чуть меньше при квадратичной классифицирующей функции
Попарный дискриминантный анализ.
partimat(new_quality ~ ., df2, method="qda")
Классификация. Априорные вероятности равны
Сравним между собой только QDA и LDA. Будем использовать априорные вероятности, равные между собой.
LDA
Y2 <- df1[-c(2,7,8)]
wine.lda <- lda(Y2, quality2, prior = c(1/3,1/3,1/3))
wine.lda## Call:
## lda(Y2, quality2, prior = c(1/3, 1/3, 1/3))
##
## Prior probabilities of groups:
## 1 2 3
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## volatile.acidity total.sulfur.dioxide pH sulphates alcohol
## 1 -0.3824288 3.262241 3.384127 -0.5654098 2.320054
## 2 -0.6701275 3.666288 3.311296 -0.4613595 2.323288
## 3 -0.9628430 3.300454 3.288802 -0.3125171 2.440156
##
## Coefficients of linear discriminants:
## LD1 LD2
## volatile.acidity -1.6472502 1.2650775
## total.sulfur.dioxide 0.1528035 -1.0232879
## pH -1.6301057 0.8590087
## sulphates 1.4154567 0.3363287
## alcohol 6.8009194 5.2880986
##
## Proportion of trace:
## LD1 LD2
## 0.8836 0.1164wine.ldap <- predict(wine.lda, Y2)$class
t1 <- table(wine.ldap, quality2)
t1## quality2
## wine.ldap 1 2 3
## 1 43 340 20
## 2 15 693 21
## 3 5 286 176Доля ошибок
1-signif(diag(prop.table(t1, 2)),3)## 1 2 3
## 0.317 0.475 0.189Ошибок в классификации 1 и 3 групп стало значительно меньше (в 3 раза), но при этом существенно увеличилось количество ошибок в классификации 2 группы (в 5 раз).
QDA
wine.qda <- qda(Y2, quality2, prior = c(1/3,1/3,1/3))
wine.qda## Call:
## qda(Y2, quality2, prior = c(1/3, 1/3, 1/3))
##
## Prior probabilities of groups:
## 1 2 3
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## volatile.acidity total.sulfur.dioxide pH sulphates alcohol
## 1 -0.3824288 3.262241 3.384127 -0.5654098 2.320054
## 2 -0.6701275 3.666288 3.311296 -0.4613595 2.323288
## 3 -0.9628430 3.300454 3.288802 -0.3125171 2.440156wine.qdap <- predict(wine.qda, Y2)$class
t2 <- table(wine.qdap, quality2)
t2## quality2
## wine.qdap 1 2 3
## 1 43 286 16
## 2 18 734 25
## 3 2 299 176Доля ошибок
1-signif(diag(prop.table(t2, 2)),3)## 1 2 3
## 0.317 0.444 0.189Качество классификации при квадратичной функции не изменилось по сравнению с линейной функцией