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

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

Необходимо построить модель на основе 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. В текст отчёта включить постановку задачи и ответы на вопросы задания.

Вариант - 13 (Упражнение 3)

13 678 Glass{mlbench} – химический состав разных типов стекла Type 2
(1 – наличие признака, все остальные – отсутствие) все остальные Логистическая регрессия, LDA

Пакеты:

library('e1071')     # SVM
library('ROCR')      # ROC-кривые
library('mlbench')   # данные Glass
library('ISLR')
library('GGally')
library('MASS')
data(Glass)
head(Glass)
##        RI    Na   Mg   Al    Si    K   Ca Ba   Fe Type
## 1 1.52101 13.64 4.49 1.10 71.78 0.06 8.75  0 0.00    1
## 2 1.51761 13.89 3.60 1.36 72.73 0.48 7.83  0 0.00    1
## 3 1.51618 13.53 3.55 1.54 72.99 0.39 7.78  0 0.00    1
## 4 1.51766 13.21 3.69 1.29 72.61 0.57 8.22  0 0.00    1
## 5 1.51742 13.27 3.62 1.24 73.08 0.55 8.07  0 0.00    1
## 6 1.51596 12.79 3.61 1.62 72.97 0.64 8.07  0 0.26    1

Зададим ядро генератора случайных чисел и объём обучающей выборки.

my.seed <- 678        # ядро генерации
train.percent <- 0.75 # доля обучающей выборки

Исходные данные: набор Glass (химический состав разных типов стекла).

Type1 <- rep(0, length(Glass$Type)) # создание вектора Type1
Glass <- cbind(Glass, Type1)        # присоединение Type1 к фрейму Glass

# замена в переменной Type: если Type = 2 означает наличие признака (1), остальные - отсутствие(0)
for(i in 1:length(Glass$Type)) {if (Glass$Type[i] == 2) {Glass$Type1[i] = 1}}

# определение долей
table(Glass$Type1) / sum(table(Glass$Type1))
## 
##         0         1 
## 0.6448598 0.3551402

Доля наименьшего класса, в данном случае 0.355, это ошибка нулевого классификатора: если бы мы прогнозировали Type = 2 для всех наблюдений, ровно в такой доле случаев мы бы ошиблись. Точность моделей целесообразно будет сравнивать с этой величиной.

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

# отбираем наблюдения в обучающую выборку 
set.seed(my.seed)
inTrain <- sample(seq_along(Glass$Type1),
                  nrow(Glass)*train.percent)
xtrain <- Glass[inTrain, c(-10, -11)]
xtest <- Glass[-inTrain, c(-10, -11)]
ytrain <- Glass[inTrain, 11]
ytest <- Glass[-inTrain, 11]

#обучающая выборка
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 cost = 3.5
## kernel = linear        0.759      0.759    0.759      0.759    0.759      0.759
## kernel = polynomial    0.741      0.741    0.722      0.704    0.704      0.704
##                     cost = 4 cost = 4.5 cost = 5 cost = 5.5 cost = 6 cost = 6.5
## kernel = linear        0.759      0.759    0.759      0.759    0.759      0.759
## kernel = polynomial    0.704      0.722    0.722      0.722    0.722      0.704
##                     cost = 7 cost = 7.5 cost = 8 cost = 8.5 cost = 9 cost = 9.5
## kernel = linear        0.759      0.778    0.815      0.815    0.815      0.815
## kernel = polynomial    0.704      0.704    0.685      0.685    0.685      0.685
##                     cost = 10 cost = 10.5 cost = 11 cost = 11.5 cost = 12
## kernel = linear         0.815       0.815     0.815       0.815     0.815
## kernel = polynomial     0.685       0.667     0.685       0.685     0.685
##                     cost = 12.5 cost = 13 cost = 13.5 cost = 14 cost = 14.5
## kernel = linear           0.815     0.778       0.796     0.778       0.778
## kernel = polynomial       0.685     0.685       0.685     0.685       0.685
##                     cost = 15 cost = 15.5 cost = 16 cost = 16.5 cost = 17
## kernel = linear         0.778       0.778     0.778       0.778     0.778
## kernel = polynomial     0.704       0.704     0.704       0.704     0.704
##                     cost = 17.5 cost = 18 cost = 18.5 cost = 19 cost = 19.5
## kernel = linear           0.778     0.778       0.778     0.778       0.778
## kernel = polynomial       0.704     0.704       0.704     0.704       0.704
##                     cost = 20
## kernel = linear         0.778
## kernel = polynomial     0.704

Из полученных результатов видно, что оптимальной формой ядерной функции будет линейная модель.

Оптимальное значение настроечного параметра

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

# список опорных векторов
svmfit$index
##   [1]   1   2   3   4   5   7   8  11  13  17  19  20  22  30  32  34  35  38
##  [19]  40  41  44  46  48  49  52  54  57  61  62  67  69  74  75  78  84  87
##  [37]  97 100 101 104 105 108 109 112 116 117 119 120 122 124 127 129 132 134
##  [55] 138 139 140 144 145 146 149 150 152 154 155 157   9  10  14  15  16  18
##  [73]  21  24  26  27  28  29  33  36  39  45  47  50  51  55  56  58  59  60
##  [91]  66  68  70  71  72  76  77  79  80  81  82  83  86  88  89  90  91  93
## [109]  94  96  98  99 102 107 110 111 114 115 121 123 126 136 142 143 147 148
## [127] 158 159
# сводка по модели
summary(svmfit)
## 
## Call:
## svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10, cale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  10 
## 
## Number of Support Vectors:  128
## 
##  ( 66 62 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1
# делаем перекрёстную проверку, изменяя штраф (аргумент cost)
set.seed(my.seed)
tune.out <- tune(svm, y ~ ., data = dat, kernel = "linear",
                 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
##  0.001
## 
## - best performance: 0.39375 
## 
## - Detailed performance results:
##    cost   error dispersion
## 1 1e-03 0.39375  0.1414520
## 2 1e-01 0.44375  0.1426741
## 3 1e+00 0.48750  0.0922331
## 4 5e+00 0.46875  0.1112196
## 5 1e+01 0.47500  0.1185854
## 6 1e+02 0.41875  0.1064337
# по перекрестной проверке наилучший настроечный параметр cost =0.01

Лучшая модель на всей обучающей выборке

# лучшая модель -- с минимальной ошибкой
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 = "linear")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  0.001 
## 
## Number of Support Vectors:  131
## 
##  ( 68 63 )
## 
## 
## 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 97 63
##        1  0  0
# чувствительность
TPR <- round(tbl1[2,2]/sum(tbl1[2,]),3)  
TPR
## [1] NaN
# специфичность
SPC <- round(tbl1[1,1]/sum(tbl1[1,]),3)  
SPC
## [1] 0.606
# функция построения 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 = "linear",  cost = 0.01, probability = T)

# матрица неточностей на обучающей (p = 0.01)
table(Predicts = predict(svmfit.opt, dat), 
             Fact = dat$y)
##         Fact
## Predicts  0  1
##        0 97 63
##        1  0  0
# прогноз вероятностей, на основе которых присваивается класс
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 41 13
##        1  0  0
# точность
ACC <- round(sum(diag(tbl2))/sum(tbl2),3)  
ACC
## [1] 0.759
# чувствительность
TPR <- round(tbl2[2,2]/sum(tbl2[2,]),3)  
TPR
## [1] NaN
# специфичность
SPC <- round(tbl2[1,1]/sum(tbl2[1,]),3)  
SPC
## [1] 0.759
# ROC-кривая для тестовой выборки
rocplot(fitted.prob, dat.te$y, main = "Test Data")
# прямая случайного классификатора
abline(a = 0, b = 1, lty = 3, lwd = 2)

Как видно из графиков ROC-кривых, и для обучающей, и для тестовой выборок значение AUC менее 0.5, а значит классификатор действует с точностью до наоборот: если положительные классификации назвать отрицательными и наоборот, классификатор будет работать лучше. Учтём это при сравнении моделей на тестовой выборке.

Сравнение моделей (логистическая регрессия, LDA, SVM) на тестовой выборке

# логистическая регрессия
model.logit <- glm(y ~ ., data = dat, family = 'binomial')
summary(model.logit)
## 
## Call:
## glm(formula = y ~ ., family = "binomial", data = dat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7155  -0.9524  -0.5824   1.1275   1.9078  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) 434.8094   323.1606   1.345   0.1785  
## x.RI         11.7999   182.2109   0.065   0.9484  
## x.Na         -4.9099     1.9741  -2.487   0.0129 *
## x.Mg         -4.0925     2.0304  -2.016   0.0438 *
## x.Al         -3.3217     2.0991  -1.582   0.1136  
## x.Si         -4.5532     1.9893  -2.289   0.0221 *
## x.K          -4.8063     2.1051  -2.283   0.0224 *
## x.Ca         -4.2000     2.0821  -2.017   0.0437 *
## x.Ba         -5.2813     2.0892  -2.528   0.0115 *
## x.Fe         -0.6644     1.9201  -0.346   0.7293  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 214.53  on 159  degrees of freedom
## Residual deviance: 190.59  on 150  degrees of freedom
## AIC: 210.59
## 
## Number of Fisher Scoring iterations: 5
# прогноз: вероятности принадлежности классу Type = 2
p.logit <- predict(model.logit, dat.te, 
                  type = 'response')

Forecast1 <- factor(ifelse(p.logit > 0.5, 2, 1),
                  levels = c(1, 2),
                  labels = c('0', '1'))

# считаем 1-SPC и TPR для всех вариантов границы отсечения
x1 <- NULL    # для (1 - SPC)
y1 <- NULL    # для TPR

# заготовка под матрицу неточностей
tbl1 <- as.data.frame(matrix(rep(0, 4), 2, 2))
rownames(tbl1) <- c('fact.0', 'fact.1')
colnames(tbl1) <- c('predict.0', 'predict.1')

# цикл по вероятностям отсечения
for (p in seq(0, 1, length = 501)){
    # прогноз
    Forecast1 <- factor(ifelse(p.logit > p, 2, 1),
                        levels = c(1, 2),
                        labels = c('0', '1'))

    # фрейм со сравнением факта и прогноза
    df.compare <- data.frame(Fact = dat.te$y, Forecast = Forecast1)

    # заполняем матрицу неточностей
    tbl1[1, 1] <- nrow(df.compare[df.compare$Fact == '0' & df.compare$Forecast == '0', ])
    tbl1[2, 2] <- nrow(df.compare[df.compare$Fact == '1' & df.compare$Forecast == '1', ])
    tbl1[1, 2] <- nrow(df.compare[df.compare$Fact == '0' & df.compare$Forecast == '1', ])
    tbl1[2, 1] <- nrow(df.compare[df.compare$Fact == '1' & df.compare$Forecast == '0', ])

    # считаем характеристики
    TPR <- tbl1[2, 2] / sum(tbl1[2, ])
    y1 <- c(y1, TPR)
    SPC <- tbl1[1, 1] / sum(tbl1[1, ])
    x1 <- c(x1, 1 - SPC)}

# LDA
model.lda <- lda(y ~ ., data = dat)

# прогноз: вероятности принадлежности классу Type = 2
p.lda <- predict(model.lda, dat.te, type = 'response')

x2 <- NULL    # для (1 - SPC)
y2 <- NULL    # для TPR

# заготовка под матрицу неточностей
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.lda$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))

# кривая (логистическая регрессия)
plot(x1, y1, type = 'l', col = 'blue', lwd = 3,
     xlab = '(1 - SPC)', ylab = 'TPR', 
     xlim = c(0, 1), ylim = c(0, 1), main = 'Тестовая выборка')

# кривая (LDA)
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', 'SVM (обр.)'), lty = 1, col = c('blue', 'red', 'green'))

Сравнивая ROC-кривые, полученные на тестовой выборке, видно, что LDA-модель обладает большей предсказательной способностью, чем логистическая регрессия и SVM(обр.).