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

Практика 9

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

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

Пример на данных по автомобилям: Auto {ISLR}.

Построим график разброса наблюдений в пространств предикторов, класс показан цветом (‘Yes’, ‘No’).

Зависимая переменная high.mpg, объсняющие: displacement, horsepower.

Переменная high.mpg – высокое значение mpg (сколько автомобиль проходит на галлоне топлива).

Перменная high.mpg принимает значение Yes, если mpg > 23, иначе - Nо.

library('e1071')     # SVM
library('ROCR')      # ROC-кривые
library('ISLR')      # данные по экспрессии генов
library('GGally')
my.seed <- 1
attach(Auto) 
head(Auto)
str(Auto)
## 'data.frame':    392 obs. of  9 variables:
##  $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
##  $ cylinders   : num  8 8 8 8 8 8 8 8 8 8 ...
##  $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
##  $ horsepower  : num  130 165 150 150 140 198 220 215 225 190 ...
##  $ weight      : num  3504 3693 3436 3433 3449 ...
##  $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
##  $ year        : num  70 70 70 70 70 70 70 70 70 70 ...
##  $ origin      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
# новая переменная 
high.mpg <- ifelse(mpg < 23, 'No', 'Yes') 
# присоединяем к таблице данных 
Auto <- data.frame(Auto, high.mpg) 
high.mpg <- as.factor(high.mpg)
ggp <- ggpairs(Auto[, c('high.mpg', 'displacement', 'horsepower')], 
               mapping = ggplot2::aes(color = high.mpg))
print(ggp, progress = F)

# таблица с данными, отклик — фактор 
dat <- data.frame(displacement, horsepower  , high.mpg) 
# обучающая выборка 
train <- sample(1:nrow(dat), nrow(dat)/2)
# SVM с радиальным ядром и маленьким cost
svmfit <- svm(high.mpg ~ ., data = dat[train, ], kernel = "radial", 
              gamma = 1, cost = 1)
plot(svmfit, dat[train, ])

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

summary(svmfit)
## 
## Call:
## svm(formula = high.mpg ~ ., data = dat[train, ], kernel = "radial", 
##     gamma = 1, cost = 1e+05)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1e+05 
##       gamma:  1 
## 
## Number of Support Vectors:  50
## 
##  ( 28 22 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  No Yes
# перекрёстная проверка
set.seed(my.seed)
tune.out <- tune(svm, high.mpg ~ ., 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     4
## 
## - best performance: 0.09657895 
## 
## - Detailed performance results:
##     cost gamma      error dispersion
## 1  1e-01   0.5 0.10710526 0.05649627
## 2  1e+00   0.5 0.10184211 0.05361723
## 3  1e+01   0.5 0.10710526 0.06170415
## 4  1e+02   0.5 0.11210526 0.05850555
## 5  1e+03   0.5 0.11736842 0.04887469
## 6  1e-01   1.0 0.11210526 0.06259742
## 7  1e+00   1.0 0.10184211 0.05361723
## 8  1e+01   1.0 0.10710526 0.06170415
## 9  1e+02   1.0 0.11210526 0.05850555
## 10 1e+03   1.0 0.11236842 0.05768236
## 11 1e-01   2.0 0.11236842 0.06374061
## 12 1e+00   2.0 0.10710526 0.06170415
## 13 1e+01   2.0 0.10157895 0.04720896
## 14 1e+02   2.0 0.11763158 0.05390349
## 15 1e+03   2.0 0.10710526 0.05542374
## 16 1e-01   3.0 0.12236842 0.07309971
## 17 1e+00   3.0 0.10184211 0.06407773
## 18 1e+01   3.0 0.11210526 0.05850555
## 19 1e+02   3.0 0.11236842 0.05207380
## 20 1e+03   3.0 0.13842105 0.06502125
## 21 1e-01   4.0 0.12763158 0.08155772
## 22 1e+00   4.0 0.09657895 0.05578354
## 23 1e+01   4.0 0.12789474 0.05668866
## 24 1e+02   4.0 0.11763158 0.05390349
## 25 1e+03   4.0 0.14342105 0.06364638

Построим матрицу неточностей для прогноза по лучшей модели и рассчитаем MSE.

# матрица неточностей для прогноза по лучшей модели
matrix <- table(true = dat[-train, "high.mpg"], 
      pred = predict(tune.out$best.model, newdata = dat[-train, ]))
bestmod <- tune.out$best.model
summary(bestmod)
## 
## Call:
## best.tune(method = svm, train.x = high.mpg ~ ., data = dat[train, 
##     ], ranges = list(cost = c(0.1, 1, 10, 100, 1000), gamma = c(0.5, 
##     1, 2, 3, 4)), kernel = "radial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  4 
## 
## Number of Support Vectors:  75
## 
##  ( 35 40 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  No Yes
#MSE
sum(diag(matrix))/sum(matrix) 
## [1] 0.9030612

ROC-кривые

# функция построения ROC-кривой: pred -- прогноз, truth -- факт
rocplot <- function(pred, truth, ...){
    predob = prediction(pred, truth)
    perf = performance(predob, "tpr", "fpr")
    plot(perf,...)}
# последняя оптимальная модель
svmfit.opt <- svm(high.mpg ~ ., data = dat[train, ], 
                  kernel = "radial", gamma = 2, cost = 0.1, decision.values = T)
# количественные модельные значения, на основе которых присваивается класс
fitted <- attributes(predict(svmfit.opt, dat[train, ],
                             decision.values = TRUE))$decision.values
# график для обучающей выборки
par(mfrow = c(1, 2))
rocplot(fitted, dat[train, "high.mpg"], main = "Training Data")
# более гибкая модель (gamma выше)
svmfit.flex = svm(high.mpg ~ ., data = dat[train, ], kernel = "radial", 
                  gamma = 50, cost = 1, decision.values = T)
fitted <- attributes(predict(svmfit.flex, dat[train, ], 
                             decision.values = T))$decision.values
rocplot(fitted, dat[train,"high.mpg"], add = T, col = "red")
# график для тестовой выборки
fitted <- attributes(predict(svmfit.opt, dat[-train, ], 
                             decision.values = T))$decision.values
rocplot(fitted, dat[-train, "high.mpg"], main = "Test Data")
fitted <- attributes(predict(svmfit.flex, dat[-train, ], 
                             decision.values = T))$decision.values
rocplot(fitted, dat[-train, "high.mpg"], add = T, col = "red")

Построенные ROC-кривые показывают большое количество точных предсказаний, что говорит о правильности выбора объясняющих перменных, однако более гибкая модель (красный цвет) хорошо работает с обучающей выборкой, но качество модели немного падает на тестовых данных, в случае другой модели качество почти не изменяется. За наилучшую и примем её- последнюю оптимальную модель.