Данные

Цветы Ирисы трех сортов.

df <- iris
df[,5] <- as.factor(df[,5]) #Изменяем тип 5-ого признака с численного на факторный, чтобы было легче дальше с ним работать.
#df <- df[,-1] #Удаляем первый столбец отвечающий за номер индивида.
head(df)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
library(matrixStats) #Позволяет делать разные прикольные штуки с матрицами - например, считать максимум по столбцам.
df.active <- df[, 1:4]
apply(df.active, 2, function(x) sum(x > 0))
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##          150          150          150          150
colMaxs(as.matrix(df.active))
## [1] 7.9 4.4 6.9 2.5

Классификация

1) Train and test

library("caTools") #Библиотека, позволяющая разбивать данные на тренировочные и тестовые случайным образом.
set.seed(13)

#dfo <- df[-c(10,77,132),] #Заранее убираем аутлаеры.
dfo <- df

sample <- sample.split(dfo, SplitRatio = 0.75) 
train <- subset(dfo, sample == TRUE)
test <- subset(dfo, sample == FALSE)
#Разбиваем данные на тренировочную и тестовую часть в отношении 3 к 1.

2) Pairs, outliers and PCA

panel.hist <- function(x, ...) 
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
pairs(train[,c(1:5)], diag.panel = panel.hist)

#Хотим посмотреть на данные
train.active <- train[,1:4]
library("FactoMineR") 
library("factoextra")
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
res2 <- PCA(train.active, scale.unit = TRUE, graph = FALSE)

fviz_pca_var(res2)

fviz_pca_biplot(res2,
                col.ind = train$Species,
                gradient.cols = c("red", "blue"),
                label = "var",
                col.var = "black",
                legend.title = "Results")

MANOVA

train.manova <- manova(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, data = train) #Смотрим на значимость дискриминации, то есть на то, отличаются ли средние по признакам между разными группами.

summary(train.manova, 'Wilks')
##           Df    Wilks approx F num Df den Df    Pr(>F)    
## Species    2 0.023954   114.69      8    168 < 2.2e-16 ***
## Residuals 87                                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(train.manova, 'Roy')
##           Df    Roy approx F num Df den Df    Pr(>F)    
## Species    2 30.062   638.83      4     85 < 2.2e-16 ***
## Residuals 87                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Выше показано, как получать результаты по разным критериям. В данном случае, значения p-value будут одинаковы для всех трех, так как данные критерии эквивалентны при сравнении двух групп.

LDA + error tables

library("MASS")

Проводим классификацию по всем данным.

full.lda <- lda(dfo[,1:4], dfo[,5]) #строим LDA
full.ldap <- predict(full.lda, dfo[,1:4]) #предсказываем по полученной классификации
#Далее строим таблицу ошибок классификации
full.ldapc <- full.ldap$class 
ct <- table(full.ldapc, dfo[,5])
ct
##             
## full.ldapc   setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         1
##   virginica       0          2        49
diag(prop.table(ct, 2)) #Столбцы - настоящие классы, строки - предсказанные классы.
##     setosa versicolor  virginica 
##       1.00       0.96       0.98
#Все аналогично, но теперь смотрим на таблицу ошибок классификации по train data на train и test отдельно.
train.lda <- lda(train[,1:4], train[,5]) 
train.ldap <- predict(train.lda, test[,1:4])
train.ldapc <- train.ldap$class
ct <- table(train.ldapc, test[,5])
ct
##             
## train.ldapc  setosa versicolor virginica
##   setosa         20          0         0
##   versicolor      0         19         1
##   virginica       0          1        19
diag(prop.table(ct, 2))
##     setosa versicolor  virginica 
##       1.00       0.95       0.95
train.ldap <- predict(train.lda, train[,1:4])
train.ldapc <- train.ldap$class
ct <- table(train.ldapc, train[,5])
ct
##             
## train.ldapc  setosa versicolor virginica
##   setosa         30          0         0
##   versicolor      0         29         1
##   virginica       0          1        29
diag(prop.table(ct, 2))
##     setosa versicolor  virginica 
##  1.0000000  0.9666667  0.9666667
#Так как данных мало, то лучше не делить на части, а использовать кросс-валидацию.
cv.lda <- lda(dfo[,1:4], dfo[,5], CV = T)
cv.ldac <- cv.lda$class
ct <- table(cv.ldac, dfo[,5])
ct
##             
## cv.ldac      setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         1
##   virginica       0          2        49
diag(prop.table(ct, 2))
##     setosa versicolor  virginica 
##       1.00       0.96       0.98

Canonical Analysis

#Этот блок будет работать, если число групп больше 2 (надо бы стрелки надписать именами признаков, иначе непонятно, какой от них толк)
#Это та плоскость, на которой различие максимально видно.
library(candisc)
## Loading required package: car
## Warning: package 'car' was built under R version 4.0.3
## Loading required package: carData
## Loading required package: heplots
## 
## Attaching package: 'candisc'
## The following object is masked from 'package:stats':
## 
##     cancor
result.candisc <- candisc(train.manova)
ggplot(result.candisc$scores) +
    geom_point(aes(x = Can1, y = Can2, color = Species)) +
    stat_ellipse(aes(x = Can1, y = Can2, color = Species)) +
    geom_segment(as.data.frame(result.candisc$structure), mapping = aes(x = 0, y = 0, xend = 6.5 * Can1, yend = 6.5 * Can2), arrow = arrow(angle = 10, type = "closed"), size = 0.5)

Exploratory graph

#Две следующие функции должны строить двумерные графики классификации, но у меня библиотека почему-то не запустилась. NE: Это из-за специфики данных.
#library(klaR)
#plot(partimat(Species ~ ., dfo, method="lda")) 

#Следующий код имеет смысл использовать, только если всего две группы и предыдущий скаттерплот не построить.
ldahist(train.ldap$x[,1], g = test[,5]) #Показывает гистограмму значений индивидов по первой канонической переменной. Здесь как раз видно, почему четыре индивида были неправильно распределены классификацией (в группе 1 - два левых столбца = 3 индивида, в группе ноль самые правый столбец)

ldahist(full.ldap$x[,1], g = full.ldapc) #Аналогично для классификации по всем данным.

ROC

#Это работает только для двух классов. Поэтому быстро повторим анализ, убрав вид цветка SETOSA
library(FSA)
## Warning: package 'FSA' was built under R version 4.0.3
## ## FSA v0.8.31. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.
## 
## Attaching package: 'FSA'
## The following object is masked from 'package:car':
## 
##     bootCase
two.dfo <- Subset(dfo,Species!="setosa")
two.lda <- lda(two.dfo[,1:4], two.dfo[,5]) 
two.ldap <- predict(two.lda, two.dfo[,1:4])
two.cv.lda <- lda(two.dfo[,1:4], two.dfo[,5], CV = T)

library(ROCR) #Библиотека позволяющая строить ROC-кривые и считать AUC

pred <- prediction(two.ldap$posterior[,2], two.dfo[,5]) #Переводит данные в стандартизированный формат. Первый параметр - предсказания, второй - настоящая классификация.
perf <- performance(pred, "tpr", "fpr") #Строит ROC-кривую
plot(perf, colorize = TRUE) #рисует ее
AUC.ROCR <- performance(pred, "auc") #Считает AUC
abline(a = 0, b = 1)
text(x = .25, y = .65, paste("AUC = ", round(AUC.ROCR@y.values[[1]],5))) #Добавляет значение AUC на график

#Аналогично делаем для кросс-валидации

predcv <- prediction(two.cv.lda$posterior[,2], two.dfo[,5])
perfcv <- performance(predcv, "tpr", "fpr")
plot(perfcv, colorize = TRUE)
AUCcv <- performance(predcv, "auc")
abline(a = 0, b = 1)
text(x = .25, y = .65, paste("AUC = ", round(AUCcv@y.values[[1]],5)))

Рисуем ROC-кривые на одной картинке. Сверху находится ROC-кривая классификации на полных данных, что проверяется по AUC

plot(perfcv, colorize = TRUE)
par(new = T)
plot(perf, colorize = TRUE)
text(x = .25, y = .55, paste("AUC_cv = ", round(AUCcv@y.values[[1]],5)))
text(x = .25, y = .65, paste("AUC_full = ", round(AUC.ROCR@y.values[[1]],4)))

Аналогично, ROC можно построить на тестовых данных, но из-за малого количества наблюдений, на нее можно не смотреть.

QDA

#Все делается аналогичным образом, только буква l в lda заменяется на q.
full.qda <- qda(dfo[,1:4], dfo[,5])
full.qdap <- predict(full.qda, dfo[,1:4])
full.qdapc <- full.qdap$class
ct <- table(full.qdapc, dfo[,5])
ct
##             
## full.qdapc   setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         1
##   virginica       0          2        49
diag(prop.table(ct, 2))
##     setosa versicolor  virginica 
##       1.00       0.96       0.98

Далее мы не будем повторять код, но нужно посмотреть на результат, построить ROC-кривую, сравнить результаты LDA и QDA