Математическое моделирование

Практика 9

Машины опорных векторов

В практических примерах ниже показано как:

  • обучить модель классификатора на опорных векторах
  • обучить модель SVM с различными формами ядерной функции
  • строить ROC-кривые встроенными функциями R

Модели: SVM
Данные: Khan {ISLR}

Подробные комментарии к коду лабораторных см. в [1], глава 9.

library('e1071')     # SVM
library('ROCR')      # ROC-кривые
library('ISLR')      # данные по экспрессии генов

my.seed <- 1

Классификатор на опорных векторах

Cгенерированные данные: два линейно неразделимых класса

# создаём наблюдения
set.seed(my.seed)
x <- matrix(rnorm(20*2), ncol = 2)
y <- c(rep(-1, 10), rep(1, 10))
x[y == 1, ] <- x[y == 1, ] + 1

# данные не разделяются линейно
plot(x, pch = 19, col = (3 - y)) 

# таблица с данными, отклик -- фактор
dat <- data.frame(x = x, y = as.factor(y))

# классификатор на опорных векторах с линейной границей
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 10, scale = FALSE)

# на графике опорные наблюдения показаны крестиками
plot(svmfit, dat)

# список опорных векторов
svmfit$index
## [1]  1  2  5  7 14 16 17
# сводка по модели
summary(svmfit)
## 
## Call:
## svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10, scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  10 
## 
## Number of Support Vectors:  7
## 
##  ( 4 3 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  -1 1
# уменьшаем штрафной параметр
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 0.1, scale = FALSE)
plot(svmfit, dat)

svmfit$index
##  [1]  1  2  3  4  5  7  9 10 12 13 14 15 16 17 18 20
# делаем перекрёстную проверку, изменяя штраф (аргумент cost)
set.seed(my.seed)
tune.out <- tune(svm, y ~ ., data = dat, kernel = "linear",
                 ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##   0.1
## 
## - best performance: 0.05 
## 
## - Detailed performance results:
##    cost error dispersion
## 1 1e-03  0.55  0.4377975
## 2 1e-02  0.55  0.4377975
## 3 1e-01  0.05  0.1581139
## 4 1e+00  0.15  0.2415229
## 5 5e+00  0.15  0.2415229
## 6 1e+01  0.15  0.2415229
## 7 1e+02  0.15  0.2415229
# лучшая модель -- с минимальной ошибкой
bestmod <- tune.out$best.model
summary(bestmod)
## 
## Call:
## best.tune(method = svm, train.x = y ~ ., data = dat, ranges = list(cost = c(0.001, 
##     0.01, 0.1, 1, 5, 10, 100)), kernel = "linear")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  0.1 
## 
## Number of Support Vectors:  16
## 
##  ( 8 8 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  -1 1
# генерируем контрольные данные
xtest <- matrix(rnorm(20*2), ncol = 2)
ytest <- sample(c(-1,1), 20, rep = TRUE)
xtest[ytest == 1, ] <- xtest[ytest == 1, ] + 1
testdat <- data.frame(x = xtest, y = as.factor(ytest))

# делаем прогноз по лучшей модели
ypred <- predict(bestmod, testdat)

# матрица неточностей
table(Прогноз = ypred, Факт = testdat$y)
##        Факт
## Прогноз -1 1
##      -1  9 1
##      1   2 8
# прогноз по модели с cost = 0.01
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = .01, scale = FALSE)
ypred <- predict(svmfit, testdat)
# матрица неточностей
table(Прогноз = ypred, Факт = testdat$y)
##        Факт
## Прогноз -1  1
##      -1 11  6
##      1   0  3

Сгенерированные данные: два линейно разделимых класса

# создаём наблюдения
x[y == 1, ] <- x[y == 1, ] + 0.5
plot(x, col = (y+5)/2, pch = 19)

# таблица с данными, отклик -- фактор
dat <- data.frame(x = x, y = as.factor(y))

# очень большой cost (маленький зазор, высокая точность классификации)
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 1e5)
summary(svmfit)
## 
## Call:
## svm(formula = y ~ ., data = dat, kernel = "linear", cost = 1e+05)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1e+05 
## 
## Number of Support Vectors:  3
## 
##  ( 1 2 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  -1 1
plot(svmfit, dat)

svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 1)
summary(svmfit)
## 
## Call:
## svm(formula = y ~ ., data = dat, kernel = "linear", cost = 1)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1 
## 
## Number of Support Vectors:  7
## 
##  ( 4 3 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  -1 1
plot(svmfit,dat)

Машина опорных векторов

Сгенерированные данные: нелинейная граница между классами

# создаём наблюдения
set.seed(my.seed)
x <- matrix(rnorm(200*2), ncol = 2)
x[1:100, ] <- x[1:100, ] + 2
x[101:150, ] <- x[101:150, ] - 2
y <- c(rep(1, 150), rep(2, 50))

# таблица с данными, отклик -- фактор 
dat <- data.frame(x = x, y = as.factor(y))
plot(x, col = y, pch = 19)

# обучающая выборка
train <- sample(200, 100)

# SVM с радиальным ядром и маленьким cost
svmfit <- svm(y ~ ., data = dat[train, ], kernel = "radial", 
              gamma = 1, cost = 1)
plot(svmfit, dat[train, ])

summary(svmfit)
## 
## Call:
## svm(formula = y ~ ., data = dat[train, ], kernel = "radial", gamma = 1, 
##     cost = 1)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
## 
## Number of Support Vectors:  31
## 
##  ( 16 15 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  1 2
# SVM с радиальным ядром и большим cost
svmfit <- svm(y ~ ., data = dat[train, ], kernel = "radial", 
              gamma = 1, cost = 1e5)
plot(svmfit, dat[train, ])

# перекрёстная проверка
set.seed(my.seed)
tune.out <- tune(svm, y ~ ., data = dat[train, ], kernel = "radial", 
                 ranges = list(cost = c(0.1, 1, 10, 100, 1000),
                               gamma = c(0.5, 1, 2, 3, 4)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost gamma
##     1   0.5
## 
## - best performance: 0.07 
## 
## - Detailed performance results:
##     cost gamma error dispersion
## 1  1e-01   0.5  0.26 0.15776213
## 2  1e+00   0.5  0.07 0.08232726
## 3  1e+01   0.5  0.07 0.08232726
## 4  1e+02   0.5  0.14 0.15055453
## 5  1e+03   0.5  0.11 0.07378648
## 6  1e-01   1.0  0.22 0.16193277
## 7  1e+00   1.0  0.07 0.08232726
## 8  1e+01   1.0  0.09 0.07378648
## 9  1e+02   1.0  0.12 0.12292726
## 10 1e+03   1.0  0.11 0.11005049
## 11 1e-01   2.0  0.27 0.15670212
## 12 1e+00   2.0  0.07 0.08232726
## 13 1e+01   2.0  0.11 0.07378648
## 14 1e+02   2.0  0.12 0.13165612
## 15 1e+03   2.0  0.16 0.13498971
## 16 1e-01   3.0  0.27 0.15670212
## 17 1e+00   3.0  0.07 0.08232726
## 18 1e+01   3.0  0.08 0.07888106
## 19 1e+02   3.0  0.13 0.14181365
## 20 1e+03   3.0  0.15 0.13540064
## 21 1e-01   4.0  0.27 0.15670212
## 22 1e+00   4.0  0.07 0.08232726
## 23 1e+01   4.0  0.09 0.07378648
## 24 1e+02   4.0  0.13 0.14181365
## 25 1e+03   4.0  0.15 0.13540064
# матрица неточностей для прогноза по лучшей модели
table(Прогноз = predict(tune.out$best.model, newdata = dat[-train, ]),
             Факт = dat[-train, "y"])
##        Факт
## Прогноз  1  2
##       1 67  2
##       2 10 21

ROC-кривые

# функция построения ROC-кривой: pred -- прогноз, truth -- факт
rocplot <- function(pred, truth, ...){
    predob = prediction(pred, truth)
    perf = performance(predob, "tpr", "fpr")
    plot(perf,...)}

# последняя оптимальная модель
svmfit.opt <- svm(y ~ ., data = dat[train, ], 
                  kernel = "radial", gamma = 2, cost = 1, probability = T)

# матрица неточностей на обучающей (p = 0.5)
table(Прогноз = predict(svmfit.opt, dat[train, ]), 
             Факт = dat[train, ]$y)
##        Факт
## Прогноз  1  2
##       1 69  3
##       2  4 24
# прогноз вероятностей, на основе которых присваивается класс
fitted.prob <- predict(svmfit.opt, dat[train, ], type = "prob",
                       probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]

# график для обучающей выборки
par(mfrow = c(1, 2))
# ROC-кривая для первой модели
rocplot(fitted.prob, dat[train, "y"], main = "Training Data")

predob = prediction(fitted.prob, dat[train, "y"])
perf = performance(predob, "tpr", "fpr")

# более гибкая модель (gamma выше)
svmfit.flex = svm(y ~ ., data = dat[train, ], 
                  kernel = "radial", gamma = 10, cost = 1, probability = T)

# матрица неточностей на обучающей (p = 0.5)
table(Прогноз = predict(svmfit.flex, dat[train, ]), 
      Факт = dat[train, ]$y)
##        Факт
## Прогноз  1  2
##       1 72  2
##       2  1 25
# прогноз вероятностей, на основе которых присваивается класс
fitted.prob <- predict(svmfit.flex, dat[train, ], type = "prob",
                       probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]
# ROC-кривая для второй модели
rocplot(fitted.prob, dat[-train, "y"], add = T, col = "red")

# график для тестовой выборки
fitted.prob <- predict(svmfit.opt, dat[-train, ], type = "prob",
                       probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]

# матрица неточностей на тестовой (p = 0.5)
table(Прогноз = predict(svmfit.opt, dat[-train, ]), 
             Факт = dat[-train, ]$y)
##        Факт
## Прогноз  1  2
##       1 66  3
##       2 11 20
rocplot(fitted.prob, dat[-train, "y"], main = "Test Data")

fitted.prob <- predict(svmfit.flex, dat[-train, ], type = "prob",
                       probability = TRUE)
fitted.prob <- attr(fitted.prob, "probabilities")[, 2]

# матрица неточностей на тестовой (p = 0.5)
table(Прогноз = predict(svmfit.flex, dat[-train, ]), 
             Факт = dat[-train, ]$y)
##        Факт
## Прогноз  1  2
##       1 67  4
##       2 10 19
rocplot(fitted.prob, dat[-train, "y"], add = T, col = "red")

SVM с несколькими классами

# генерируем данные
set.seed(my.seed)
x <- rbind(x, matrix(rnorm(50*2), ncol = 2))
y <- c(y, rep(0, 50))
x[y == 0, 2] <- x[y == 0, 2] + 2
dat <- data.frame(x = x, y = as.factor(y))

# график и модель по методу "один против одного"
par(mfrow = c(1, 1))
plot(x, col = (y + 1))

svmfit = svm(y ~ ., data = dat, kernel = "radial", cost = 10, gamma = 1)
plot(svmfit, dat)

Сеточный поиск

Анализ данных по уровню экспрессии генов

# данные по образцам тканей четырёх типов саркомы
names(Khan)
## [1] "xtrain" "xtest"  "ytrain" "ytest"
dim(Khan$xtrain)     # обучающая выборка, предикторы
## [1]   63 2308
dim(Khan$xtest)      # тестовая выборка, предикторы
## [1]   20 2308
length(Khan$ytrain)  # обучающая выборка, отклик
## [1] 63
length(Khan$ytest)   # тестовая выборка, отклик
## [1] 20
table(Khan$ytrain)
## 
##  1  2  3  4 
##  8 23 12 20
table(Khan$ytest)
## 
## 1 2 3 4 
## 3 6 6 5
dat <- data.frame(x = Khan$xtrain, y = as.factor(Khan$ytrain))

# тестовые данные
dat.te <- data.frame(x = Khan$xtest, y = as.factor(Khan$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"
AUC
##                     cost = 1 cost = 1.5 cost = 2 cost = 2.5 cost = 3 cost = 3.5
## kernel = linear          0.9       0.90      0.9       0.90     0.90       0.90
## kernel = polynomial      0.3       0.35      0.5       0.55     0.55       0.55
##                     cost = 4 cost = 4.5 cost = 5 cost = 5.5 cost = 6 cost = 6.5
## kernel = linear         0.90       0.90     0.90       0.90     0.90       0.90
## kernel = polynomial     0.55       0.55     0.55       0.55     0.55       0.55
##                     cost = 7 cost = 7.5 cost = 8 cost = 8.5 cost = 9 cost = 9.5
## kernel = linear         0.90       0.90     0.90       0.90     0.90       0.90
## kernel = polynomial     0.55       0.55     0.55       0.55     0.55       0.55
##                     cost = 10 cost = 10.5 cost = 11 cost = 11.5 cost = 12
## kernel = linear          0.90        0.90      0.90        0.90      0.90
## kernel = polynomial      0.55        0.55      0.55        0.55      0.55
##                     cost = 12.5 cost = 13 cost = 13.5 cost = 14 cost = 14.5
## kernel = linear            0.90      0.90        0.90      0.90        0.90
## kernel = polynomial        0.55      0.55        0.55      0.55        0.55
##                     cost = 15 cost = 15.5 cost = 16 cost = 16.5 cost = 17
## kernel = linear          0.90        0.90      0.90        0.90      0.90
## kernel = polynomial      0.55        0.55      0.55        0.55      0.55
##                     cost = 17.5 cost = 18 cost = 18.5 cost = 19 cost = 19.5
## kernel = linear            0.90      0.90        0.90      0.90        0.90
## kernel = polynomial        0.55      0.55        0.55      0.55        0.55
##                     cost = 20
## kernel = linear          0.90
## kernel = polynomial      0.55

Упражнение 9

Необходимо построить модель на основе SVM для указанной в варианте зависимой переменной.

Данные взять из упражнения №3

Для модели:

1 Отложить 25% наблюдений в тестовую выборку (ядро генератора случайных чисел указано в варианте к упражнению №3).
2 На обучающей выборке (оставшихся 75% наблюдений) сравнить несколько видов ядер SVM по точности модели (AUC) методом сеточного поиска.
3 Для оптимальной формы ядерной функции на обучающей выборке подобрать оптимальное значение настроечных параметров по минимальной ошибке с перекрёстной проверкой (функция tune).
4 Подогнать лучшую модель на всей обучающей выборке. Построить ROC-кривую и рассчитать матрицу неточностей, чувствительность и специфичность.
5 Сделать прогноз по лучшей модели на тестовую выборку, оценить его качество точность по матрице неточностей, чувствительность и специфичность, построить ROC-кривую.
6 Сравнить результаты, которые дал SVM, с результатами, полученными в упражнении 3. Какой из методов оказался лучше?

Как сдавать: прислать на почту преподавателя ссылки: * на html-отчёт с видимыми блоками кода (блоки кода с параметром echo = T), размещённый на rpubs.com.
* на код, генерирующий отчёт, в репозитории на github.com. В текст отчёта включить постановку задачи и ответы на вопросы задания.

Источники

  1. Джеймс Г., Уиттон Д., Хасти Т., Тибширани Р. Введение в статистическое обучение с примерами на языке R / пер. с англ. С.Э. Мастицкого. – М.: ДМК Пресс, 2016 – 450 с. Репозиторий с примерами к книге на русском языке: https://github.com/ranalytics/islr-ru
  2. Х.Бринк, Д.Ричардс, М.Феверолф Машинное обучение. – М.: СПб.: Питер, 2018 – 336 с.