Лабораторная №9
Косован Михаил
13 05 2020
Необходимо построить модель на основе SVM для указанной в варианте зависимой переменной.
Данные взять из упражнения №3.
Для модели:
Отложить 25% наблюдений в тестовую выборку (ядро генератора случайных чисел указано в варианте к упражнению №3).
На обучающей выборке (оставшихся 75% наблюдений) сравнить несколько видов ядер SVM по точности модели (AUC) методом сеточного поиска.
Для оптимальной формы ядерной функции на обучающей выборке подобрать оптимальное значение настроечных параметров по минимальной ошибке с перекрёстной проверкой (функция tune).
Подогнать лучшую модель на всей обучающей выборке. Построить ROC-кривую и рассчитать матрицу неточностей, чувствительность и специфичность.
Сделать прогноз по лучшей модели на тестовую выборку, оценить его качество точность по матрице неточностей, чувствительность и специфичность, построить ROC-кривую.
Сравнить результаты, которые дал SVM, с результатами, полученными в упражнении 3. Какой из методов оказался лучше?
Вариант 10
Ядро для set.seed() - 345.
Данные: titanic_train{titanic} - выжившие в катастрофе Титаника.
Зависимая переменная: survived.
Объясняющие переменные: Все остальные, кроме name.
Методы: LDA, QDA.
Пакеты:
library('e1071') # SVM
library('ROCR') # ROC-кривые
library('ISLR')
library('GGally')
library('MASS')
library('titanic')
data(titanic_train)
head(titanic_train)Зададим ядро генератора случайных чисел и объем обучающей выборки.
# Зададим ядро генератора случайных чисел и объем обучающей выборки
my.seed <- 345
train.percent <- 0.75
options("ggmatrix.progress.bar" = FALSE)Исходные данные: набор titanic_train (выжившие в катастрофе Титаника)
ggp <- ggpairs(titanic_train[c(-4,-9, -11)])
print(ggp, progress = FALSE)
titanic_train <- titanic_train[-4]
table(titanic_train$Survived) / sum(table(titanic_train$Survived))##
## 0 1
## 0.6161616 0.3838384Для наименьшего класса, в данном случае 0.384, это ошибка нулевого классификатора: если бы мы прогнозировали Survived = 1 для всех наблюдений, ровно в такой доле случаев мы бы ошиблись. Точность моделей целесообразно будет сравнивать с этой величиной.
Сеточный поиск
# Убираем пропуски
titanic_train <- na.omit(titanic_train)
# Отбираем наблюдения в обучающую выборку
set.seed(my.seed)
inTrain <- sample(seq_along(titanic_train$Survived),
nrow(titanic_train)*train.percent)
xtrain <- titanic_train[inTrain, c(-2, -8, -10, -11)]
xtest <- titanic_train[-inTrain, c(-2, -8, -10, -11)]
ytrain <- titanic_train[inTrain, 2]
ytest <- titanic_train[-inTrain, 2]
# Обучающая выборка
dat <- data.frame(x = xtrain, y = as.factor(ytrain))
# Тестовые данные
dat.te <- data.frame(x = xtest, y = as.factor(ytest))
# Параметры алгоритма
kernel.grid <- c('linear', 'polynomial')
cost.grid <- seq(1, 20, by = 0.5)
AUC <- matrix(0, length(kernel.grid), length(cost.grid))
colnames(AUC) <- paste0('cost = ', cost.grid)
rownames(AUC) <- paste0('kernel = ', kernel.grid)
# SVM
for (i in 1:length(kernel.grid)) {
print(paste0('Starting ', kernel.grid[i], ' kernel'))
for (j in 1:length(cost.grid)) {
out <- svm(y ~ ., data = dat, kernel = kernel.grid[i],
cost = cost.grid[j])
# Прогноз на тестовой выборке
pred.te <- predict(out, newdata = dat.te)
# Матрица неточностей
tbl <- table(pred.te, dat.te$y)
AUC[i, j] <- sum(diag(tbl)) / sum(tbl)
}
}## [1] "Starting linear kernel"
## [1] "Starting polynomial kernel"round(AUC, 3)## cost = 1 cost = 1.5 cost = 2 cost = 2.5 cost = 3
## kernel = linear 0.804 0.804 0.804 0.804 0.804
## kernel = polynomial 0.771 0.788 0.793 0.799 0.810
## cost = 3.5 cost = 4 cost = 4.5 cost = 5 cost = 5.5
## kernel = linear 0.804 0.804 0.804 0.804 0.804
## kernel = polynomial 0.821 0.816 0.827 0.832 0.832
## cost = 6 cost = 6.5 cost = 7 cost = 7.5 cost = 8
## kernel = linear 0.804 0.804 0.804 0.804 0.804
## kernel = polynomial 0.832 0.832 0.827 0.821 0.821
## cost = 8.5 cost = 9 cost = 9.5 cost = 10 cost = 10.5
## kernel = linear 0.804 0.804 0.804 0.804 0.804
## kernel = polynomial 0.816 0.816 0.804 0.804 0.804
## cost = 11 cost = 11.5 cost = 12 cost = 12.5 cost = 13
## kernel = linear 0.804 0.804 0.804 0.804 0.804
## kernel = polynomial 0.799 0.799 0.799 0.793 0.799
## cost = 13.5 cost = 14 cost = 14.5 cost = 15
## kernel = linear 0.804 0.804 0.804 0.804
## kernel = polynomial 0.799 0.799 0.799 0.804
## cost = 15.5 cost = 16 cost = 16.5 cost = 17
## kernel = linear 0.804 0.804 0.804 0.804
## kernel = polynomial 0.804 0.804 0.810 0.810
## cost = 17.5 cost = 18 cost = 18.5 cost = 19
## kernel = linear 0.804 0.804 0.804 0.804
## kernel = polynomial 0.810 0.816 0.804 0.804
## cost = 19.5 cost = 20
## kernel = linear 0.804 0.804
## kernel = polynomial 0.804 0.799Из полученных результатов видно, что оптимальной формой ядерной функции будет полиномиальная модель.
Оптимальное значение настроечного параметра
# Классификатор на опорных векторах
svmfit <- svm(y ~ ., data = dat, kernel = "polynomial", cost = 6, cale = FALSE)
# Список опорных векторов
svmfit$index## [1] 1 3 7 11 13 19 29 30 32 36 39 41 44 49 59 65 71
## [18] 77 85 93 95 96 101 103 107 110 111 113 114 115 116 117 119 126
## [35] 133 134 136 138 142 143 145 161 167 180 184 186 187 189 196 197 200
## [52] 201 205 206 217 218 221 230 232 239 240 243 245 248 249 253 255 257
## [69] 260 268 278 288 292 295 300 305 307 311 314 316 317 320 327 334 351
## [86] 353 357 358 360 373 384 390 392 394 396 397 413 414 419 421 424 427
## [103] 428 429 431 434 442 443 445 446 453 456 474 475 476 485 487 492 494
## [120] 495 499 502 506 507 511 515 517 518 521 529 532 535 4 5 6 9
## [137] 10 15 17 18 33 35 42 46 47 62 72 78 80 86 87 88 89
## [154] 90 92 106 109 112 118 120 123 124 127 128 131 146 151 152 156 158
## [171] 160 164 165 168 170 171 173 179 182 185 191 198 202 209 220 222 223
## [188] 224 227 229 233 234 236 237 242 254 258 272 273 275 281 285 286 290
## [205] 291 296 297 299 302 306 309 323 324 325 338 347 352 354 366 369 371
## [222] 375 378 379 381 383 386 387 388 395 398 401 404 406 407 408 415 417
## [239] 430 433 440 441 448 450 452 454 455 457 458 459 460 466 467 468 470
## [256] 482 483 493 496 500 505 508 509 514 525 531# Сводка по модели
summary(svmfit)##
## Call:
## svm(formula = y ~ ., data = dat, kernel = "polynomial", cost = 6,
## cale = FALSE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 6
## degree: 3
## coef.0: 0
##
## Number of Support Vectors: 266
##
## ( 132 134 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1# Делаем перекрёстную проверку, изменяя штраф (аргумент cost)
set.seed(my.seed)
tune.out <- tune(svm, y ~ ., data = dat, kernel = "polynomial",
ranges = list(cost = c(0.001, 0.1, 1, 5, 10, 100)))
summary(tune.out)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 10
##
## - best performance: 0.1983229
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-03 0.4040182 0.08380033
## 2 1e-01 0.3237945 0.08969253
## 3 1e+00 0.2378057 0.09159166
## 4 5e+00 0.2078616 0.08229671
## 5 1e+01 0.1983229 0.06450308
## 6 1e+02 0.2058351 0.06093292# По перекрестной проверке наилучший настроечный параметр cost = 10Лучшая модель на всей обучающей выборке
# Лучшая модель -- с минимальной ошибкой
bestmod <- tune.out$best.model
summary(bestmod)##
## Call:
## best.tune(method = svm, train.x = y ~ ., data = dat, ranges = list(cost = c(0.001,
## 0.1, 1, 5, 10, 100)), kernel = "polynomial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 10
## degree: 3
## coef.0: 0
##
## Number of Support Vectors: 254
##
## ( 127 127 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1# Делаем прогноз по лучшей модели
ypred_tr <- predict(bestmod, dat)
# Матрица неточностей
tbl1 <- table(Predicts = ypred_tr, Fact = dat$y)
tbl1## Fact
## Predicts 0 1
## 0 295 46
## 1 23 171# Чувствительность
TPR <- round(tbl1[2,2]/sum(tbl1[2,]),3)
TPR## [1] 0.881# Специфичность
SPC <- round(tbl1[1,1]/sum(tbl1[1,]),3)
SPC## [1] 0.865# Функция построения ROC-кривой: pred -- прогноз, truth -- факт
rocplot <- function(pred, truth, ...){
predob = prediction(pred, truth)
perf = performance(predob, "tpr", "fpr")
plot(perf,...)}
# Оптимальная модель
svmfit.opt <- svm(y ~ ., data = dat, kernel = "polynomial", cost = 10, probability = T)
# Матрица неточностей на обучающей (p = 0.01)
table(Predicts = predict(svmfit.opt, dat),
Fact = dat$y)## Fact
## Predicts 0 1
## 0 295 46
## 1 23 171# Прогноз вероятностей, на основе которых присваивается класс
fitted.prob <- predict(svmfit.opt, dat, type = "prob", probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]
# График для обучающей выборки
# ROC-кривая для первой модели
rocplot(fitted.prob, dat[, "y"], main = "Training Data")
# Прямая случайного классификатора
abline(a = 0, b = 1, lty = 3, lwd = 2)
Лучшая модель на тестовой выборке
# График для тестовой выборки
fitted.prob <- predict(svmfit.opt, dat.te, type = "prob",
probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]
# Матрица неточностей на тестовой (p = 0.01)
tbl2 <- table(Predicts = predict(svmfit.opt, dat.te), Fact = dat.te$y)
tbl2## Fact
## Predicts 0 1
## 0 92 21
## 1 14 52# Точность
ACC <- round(sum(diag(tbl2))/sum(tbl2),3)
ACC## [1] 0.804# Чувствительность
TPR <- round(tbl2[2,2]/sum(tbl2[2,]),3)
TPR## [1] 0.788# Специфичность
SPC <- round(tbl2[1,1]/sum(tbl2[1,]),3)
SPC## [1] 0.814# ROC-кривая для тестовой выборки
rocplot(fitted.prob, dat.te$y, main = "Test Data")
# Прямая случайного классификатора
abline(a = 0, b = 1, lty = 3, lwd = 2)
Как видно из графиков ROC-кривых, и для обучающей, и для тестовой выборок значение AUC менее 0.5, а значит классификатор действует с точностью до наоборот: если положительные классификации назвать отрицательными и наоборот, классификатор будет работать лучше. Учтём это при сравнении моделей на тестовой выборке.
# LDA
model.lda <- lda(y ~ ., data = dat)
# Прогноз: вероятности принадлежности классу Survived = 1
p.lda <- predict(model.lda, dat.te, type = 'response')
Forecast1 <- factor(ifelse(p.lda$posterior[, '1'] > 0.5, 2, 1), levels = c(1, 2), labels = c('0', '1'))
# Для (1 - SPC)
x1 <- NULL
# Для TPR
y1 <- NULL
# Заготовка под матрицу неточностей
tbl2 <- as.data.frame(matrix(rep(0, 4), 2, 2))
rownames(tbl2) <- c('fact.0', 'fact.1')
colnames(tbl2) <- c('predict.0', 'predict.1')
# Цикл по вероятностям отсечения
for (p in seq(0, 1, length = 501)){
# Прогноз
Forecast1 <- factor(ifelse(p.lda$posterior[, '1'] > p, 2, 1),
levels = c(1, 2),
labels = c('0', '1'))
# фрейм со сравнением факта и прогноза
df.compare <- data.frame(Fact = dat.te$y, Forecast = Forecast1)
# Заполняем матрицу неточностей
tbl2[1, 1] <- nrow(df.compare[df.compare$Fact == '0' & df.compare$Forecast == '0', ])
tbl2[2, 2] <- nrow(df.compare[df.compare$Fact == '1' & df.compare$Forecast == '1', ])
tbl2[1, 2] <- nrow(df.compare[df.compare$Fact == '0' & df.compare$Forecast == '1', ])
tbl2[2, 1] <- nrow(df.compare[df.compare$Fact == '1' & df.compare$Forecast == '0', ])
# Считаем характеристики
TPR <- tbl2[2, 2] / sum(tbl2[2, ])
y1 <- c(y1, TPR)
SPC <- tbl2[1, 1] / sum(tbl2[1, ])
x1 <- c(x1, 1 - SPC)
}
model.qda <- qda(y ~ ., data = dat)
# Прогноз: вероятности принадлежности классу Survived = 1
p.qda <- predict(model.qda, dat.te, type = 'response')
Forecast2 <- factor(ifelse(p.qda$posterior[, '1'] > 0.5, 2, 1), levels = c(1, 2), labels = c('0', '1'))
# Для (1 - SPC)
x2 <- NULL
# Для TPR
y2 <- NULL
# Заготовка под матрицу неточностей
tbl2 <- as.data.frame(matrix(rep(0, 4), 2, 2))
rownames(tbl2) <- c('fact.0', 'fact.1')
colnames(tbl2) <- c('predict.0', 'predict.1')
# Цикл по вероятностям отсечения
for (p in seq(0, 1, length = 501)){
# Прогноз
Forecast2 <- factor(ifelse(p.qda$posterior[, '1'] > p, 2, 1),
levels = c(1, 2),
labels = c('0', '1'))
# фрейм со сравнением факта и прогноза
df.compare <- data.frame(Fact = dat.te$y, Forecast = Forecast2)
# Заполняем матрицу неточностей
tbl2[1, 1] <- nrow(df.compare[df.compare$Fact == '0' & df.compare$Forecast == '0', ])
tbl2[2, 2] <- nrow(df.compare[df.compare$Fact == '1' & df.compare$Forecast == '1', ])
tbl2[1, 2] <- nrow(df.compare[df.compare$Fact == '0' & df.compare$Forecast == '1', ])
tbl2[2, 1] <- nrow(df.compare[df.compare$Fact == '1' & df.compare$Forecast == '0', ])
# Считаем характеристики
TPR <- tbl2[2, 2] / sum(tbl2[2, ])
y2 <- c(y2, TPR)
SPC <- tbl2[1, 1] / sum(tbl2[1, ])
x2 <- c(x2, 1 - SPC)
}
# строим ROC-кривую
par(mar = c(5, 5, 1, 1))
# кривая (LDA)
plot(x1, y1, type = 'l', col = 'blue', lwd = 3,
xlab = '(1 - SPC)', ylab = 'TPR',
xlim = c(0, 1), ylim = c(0, 1), main = 'Test sample')
# кривая (QDA)
lines(x2, y2, type = 'l', col = 'red', lwd = 3)
# кривая (SVM обр.)
rocplot(-fitted.prob, dat.te$y, add = T, col = 'green')
# прямая случайного классификатора
abline(a = 0, b = 1, lty = 3, lwd = 2)
# легенда
legend('bottomright', names <- c('LDA', 'QDA', 'SVM (reverse)'), lty = 1, col = c('blue', 'red', 'green'))
Сравнивая ROC-кривые, полученные на тестовой выборке, видно, что LDA обладает большей предсказательной способностью, чем QDA и SVM(обр.).