Tutorial 5: Máquinas de Vetores de SuporteSupport Vector Machines (SVM)

Este tutorial é baseado nos exemplos do livro An Introduction to Statistical Learning, de Gareth James, Daniela Witten, Trevor Hastie e Rob Tibshirani (https://www.statlearning.com/).

Curiosamente, o livro dá a entender que o termo SVM é utilizado exclusivamente para classificadores de vetores de suporte (support vector classifier) não lineares (ou seja, baseados em funções kernel não lineares). Esse não é um consenso universal: em geral, o termo SVM é amplamente utilizado para se referir a tanto classificadores lineares quanto não lineares, bem como regressões por vetores de suporte (support vector regression) lineares e não lineares.

Classificador de Vetores de Suporte com Kernel Linear

Linear Support Vector Classifier

A biblioteca e1071 contém implementações para diversos métodos de aprendizado estatístico. Em particular, a função svm() pode ser usada para ajustar um classificador de vetores de suporte linear quando o argumento kernel = "linear" é utilizado. O argumento cost nos permite especificar o custo de uma violação da margem. Ou seja, ele é um hiperparâmetro de regularização que controla a penalidade associada à violação da margem. Quando o argumento cost é pequeno, as margens serão largas e muitos vetores de suporte estarão na margem ou a violarão. Quando o argumento cost é grande, as margens serão estreitas e haverá poucos vetores de suporte na margem ou a violando.

Utilizaremos a função svm() para ajustar um classificador de vetores de suporte linear para um dado valor do parâmetro cost. Aqui, demonstramos o uso desta função em um exemplo bidimensional para que possamos plotar a fronteira de decisão resultante. Começamos gerando as observações, que pertencem a duas classes, e verificando se as observações dessas classes são linearmente separáveis:

set.seed(1)
x <- matrix(rnorm(20 * 2), ncol = 2)
y <- c(rep(-1, 10), rep(1, 10))
x[y == 1, ] <- x[y == 1, ] + 1
plot(x, col = (3 - y), pch = 19)

Vemos que elas não são linearmente separáveis…

Em seguida, ajustamos o classificador de vetores de suporte. Note que, para que a função svm() realize a classificação (em vez de regressão baseada em SVM, ou seja, support vector regression), devemos codificar a resposta como uma variável categórica, utilizando as.factor(y). Então, criamos um data frame com a resposta codificada como um fator:

dat <- data.frame(x = x, y = as.factor(y))
library(e1071)
svmfit <- svm(y ~ ., data = dat, kernel = "linear", 
    cost = 10, scale = FALSE)

O argumento scale = FALSE informa à função svm() para não normalizar cada variável de forma a ter média \(0\) ou desvio padrão \(1\), mas dependendo da aplicação, pode-se preferir usar scale = TRUE.

Agora podemos plotar o classificador de vetores de suporte obtido:

plot(svmfit, dat)

Note que os dois argumentos para a função plot() da SVM são o output da função svm() e os dados usados (dat). A região do espaço de características que será atribuída à classe \(-1\) é mostrada em amarelo claro, e a região que será atribuída à classe \(+1\) é mostrada em vermelho.

A fronteira de decisão entre as duas classes é linear (porque usamos o argumento kernel = "linear"), embora devido à forma como a função de plotagem é implementada nesta biblioteca, a fronteira de decisão pareça um tanto irregular no gráfico. (Note que aqui a segunda característica é plotada no eixo \(x\) e a primeira característica é plotada no eixo \(y\), em contraste com o comportamento da função plot() do R.)

Os vetores de suporte são plotados como cruzes e as observações restantes são plotadas como círculos; vemos aqui que existem \(7\) vetores de suporte. Podemos determinar suas identidades da seguinte forma:

svmfit$index
[1]  1  2  5  7 14 16 17

Podemos obter algumas informações básicas sobre o classificador de vetores de suporte ajustado usando o comando summary():

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

O output acima nos informa, por exemplo, que um kernel linear foi usado com cost = 10, e que \(7\) vetores de suporte foram selecionados, \(4\) de uma classe e \(3\) da outra.

Mas o que aconteceria se usássemos um valor menor para o parâmetro de custo?

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

Quando um valor menor para o parâmetro de custo é utilizado, obtemos um número maior de vetores de suporte, porque a margem agora é mais larga. Infelizmente, a função svm() não fornece explicitamente os coeficientes da fronteira de decisão linear obtida quando o classificador de vetores de suporte é ajustado, nem fornece a largura da margem.

Mas como selecionar o “melhor” valor para o parametro cost da nossa SVM? A biblioteca e1071 inclui uma função integrada, tune(), para realizar validação cruzada. Por padrão, tune() realiza validação cruzada de dez folds em um conjunto de modelos de interesse. Para usar essa função, passamos informações relevantes sobre o conjunto de modelos que estão sob consideração. O seguinte comando indica que queremos comparar SVMs com um kernel linear usando uma variedade de valores do parâmetro cost:

set.seed(1)
tune.out <- tune(svm, y ~ ., data = dat, kernel = "linear", 
    ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)))

Podemos acessar facilmente os erros de validação cruzada para cada um desses modelos usando o comando summary():

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

Vemos que cost = 0.1 resulta na menor taxa de erro de validação cruzada. A função tune() armazena o melhor modelo obtido, que pode ser acessado da seguinte forma:

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

A função predict() pode ser usada para prever o rótulo de classe em um conjunto de observações de teste, para qualquer valor dado do parâmetro de custo. Então, começamos gerando um conjunto de dados de teste:

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))

Agora prevemos os rótulos de classe dessas observações de teste. Aqui utilizamos o melhor modelo obtido através da validação cruzada para fazer as previsões:

ypred <- predict(bestmod, testdat)
table(predict = ypred, truth = testdat$y)
       truth
predict -1 1
     -1  9 1
     1   2 8

Assim, com este valor de cost, \(17\) das observações de teste são classificadas corretamente (\(9+8=17\)) e \(3\) incorretamente (\(2+1=3\)). E se, em vez disso, tivéssemos usado cost = 0.01?

svmfit <- svm(y ~ ., data = dat, kernel = "linear", 
    cost = .01, scale = FALSE)
ypred <- predict(svmfit, testdat)
table(predict = ypred, truth = testdat$y)
       truth
predict -1  1
     -1 11  6
     1   0  3

Neste caso, \(14\) observações são classificadas corretamente (\(11 + 3 = 14\)) e \(6\) observações são classificadas incorretamente (\(6 + 0 = 6\)). Ou seja, \(3\) a mais que no caso anterior.

Agora, considere uma situação em que as duas classes são linearmente separáveis. Então, podemos encontrar um hiperplano separador usando a função svm(). Primeiro, separamos ainda mais as duas classes em nossos dados simulados para que sejam linearmente separáveis:

x[y == 1, ] <- x[y == 1, ] + 0.5
plot(x, col = (y + 5) / 2, pch = 19)

Vemos que as observações são apenas marginalmente linearmente separáveis. Agora, ajustamos o classificador de vetores de suporte e plotamos o hiperplano resultante, utilizando um valor alto de cost para que nenhuma observação seja classificada incorretamente:

dat <- data.frame(x = x, y = as.factor(y))
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)

Nenhum erro de treinamento foi cometido e apenas \(3\) vetores de suporte foram utilizados. No entanto, podemos ver na figura que a margem é muito estreita (porque as observações que não são vetores de suporte, indicadas como círculos, estão bem próximas da fronteira de decisão). Parece provável que este modelo terá um desempenho ruim em dados de teste. Agora, tentamos um valor menor de cost:

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)

Utilizando cost = 1, classificamos incorretamente uma observação de treinamento, mas também obtemos uma margem muito mais larga, e utilizamos \(7\) vetores de suporte. Parece provável que este modelo terá um desempenho melhor em dados de teste do que o modelo com cost = 1e5.

Classificador de Vetores de Suporte com Kernel Não Linear

Nonlinear Support Vector Classifier

Para ajustar um classificador de vetores de suporte não linear, utilizamos mais uma vez a função svm(). No entanto, agora usamos um valor diferente para o parâmetro kernel. Para ajustar a SVM com kernel polinomial, usamos kernel = "polynomial", e para ajustar a SVM com kernel radial, usamos kernel = "radial". No primeiro caso, também usamos o argumento degree para especificar o grau do kernel polinomial, e no último caso, usamos gamma para especificar um valor de \(\gamma\) (gama) para o kernel de base radial (revise as equações das funções kernel se necessário).

Primeiro, geramos alguns dados difíceis de serem classificados de forma linear:

set.seed(1)
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))

A plotagem dos dados deixa claro que, de fato, precisamos de uma abordagem não linear para a classificação desses dados:

plot(x, col = y)

A seguir, divididos aleatoriamente os dados em grupos de treinamento e teste. Em seguida, ajustamos os dados de treinamento usando a função svm() com um kernel radial e \(\gamma=1\):

train <- sample(200, 100)
svmfit <- svm(y ~ ., data = dat[train, ], kernel = "radial",  
    gamma = 1, cost = 1)
plot(svmfit, dat[train, ])

O gráfico mostra que a SVM resultante possui uma fronteira decididamente não linear. A função summary() pode ser usada para obter algumas informações sobre o ajuste da SVM:

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

Podemos ver na figura acima que há um número razoável de erros de treinamento neste ajuste da SVM. Se aumentarmos o valor de cost, podemos reduzir o número de erros de treinamento. No entanto, isso tem o custo de uma fronteira de decisão mais irregular que parece estar em risco de sobreajustar os dados (overfitting).

svmfit <- svm(y ~ ., data = dat[train, ], kernel = "radial", 
    gamma = 1, cost = 1e5)
plot(svmfit, dat[train, ])

Como vimos anteriormente, podemos realizar validação cruzada usando tune() para selecionar os “melhores” valores de \(\gamma\) e cost para nossa SVM com kernel radial:

set.seed(1)
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

O output acima indica que devemos adotar cost = 1 e gamma = 0.5 em nossa SVM.

Podemos visualizar as previsões no conjunto de teste do nosso modelo aplicando a função predict(). Note que para fazer isso, subdefinimos o data frame dat usando -train:

table(
    true = dat[-train, "y"], 
    pred = predict(
      tune.out$best.model, newdata = dat[-train, ]
      )
    )
    pred
true  1  2
   1 67 10
   2  2 21

Vemos que \(12 \%\) das observações de teste são classificadas incorretamente por esta SVM. Será que uma SVM com kernel polinomial se sairia melhor?

Comparação de Classificadores de Vetores de Suporte com curvas ROC

O pacote ROCR pode ser usado para produzir curvas ROC, frequentemente usadas para avaliar o desempenho de modelos de classificação binária—como os classificadores SVM deste tutorial.

Vamos ver em detalhes como curvas ROC funcionam quando falarmos sobre validação de modelos. Por hora, vamos escrever uma função simples para plotar uma curva ROC, dado um vetor contendo uma pontuação numérica para cada observação, pred, e um vetor contendo o rótulo de classe para cada observação, truth:

library(ROCR)
rocplot <- function(pred, truth, ...) {
  predob <- prediction(pred, truth)
  perf <- performance(predob, "tpr", "fpr")
  plot(perf, ...)
}

Classificadores SVM fornecem rótulos de classe para cada observação. No entanto, também é possível obter valores ajustados para cada observação, que são as pontuações numéricas usadas para obter os rótulos de classe. Por exemplo (mas sem entrar em detalhes), no caso de um classificador de vetores de suporte linear, o valor ajustado para uma observação \(X= (X_1, X_2, \ldots, X_p)^T\) assume a forma \(\hat{\beta}_0 + \hat{\beta}_1 X_1 + \hat{\beta}_2 X_2 + \cdots + \hat{\beta}_p X_p\), enquanto que para uma SVM com um kernel não linear, a equação que fornece o valor ajustado é dada por \(f(x) = \sum_{i \in \text{Support Vectors}} \alpha_i y_i K(x_i, x) + b\).

Em essência, o sinal do valor ajustado determina em qual lado da fronteira de decisão a observação se encontra. Portanto, a relação entre o valor ajustado e a previsão de classe para uma dada observação é simples: se o valor ajustado for maior que zero, a observação é atribuída a uma classe e, se for menor que zero, é atribuída à outra.

Para obter os valores ajustados para um determinado ajuste de uma SVM, usamos decision.values = TRUE ao ajustar svm(). Dessa forma, a função predict() fornecerá esses valores:

svmfit.opt <- svm(y ~ ., data = dat[train, ], 
    kernel = "radial", gamma = 2, cost = 1, 
    decision.values = T)
fitted <- attributes(
    predict(svmfit.opt, dat[train, ], decision.values = TRUE)
  )$decision.values

Com isso, podemos produzir uma curva ROC. Note que usamos os valores ajustados negativos como pertencentes à classe \(1\) e os valores positivos à classe \(2\):

par(mfrow = c(1, 2))
rocplot(-fitted, dat[train, "y"], main = "Training Data")

A SVM parece produzir previsões precisas, mas ao aumentar \(\gamma\), podemos produzir um ajuste mais flexível e melhorar ainda mais sua precisão:

rocplot(-fitted, dat[train, "y"], main = "Training Data")
svmfit.flex <- svm(y ~ ., 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, "y"], add = T, col = "red")

No entanto, essas curvas ROC são todas baseadas nos dados de treinamento. O que mais interessa é o nível de precisão das previsões quando a SVM é aplicada aos dados de teste. Quando computamos as curvas ROC usando os dados de teste, o modelo com \(\gamma=2\) parece fornecer os resultados mais precisos:

fitted <- attributes(
    predict(svmfit.opt, dat[-train, ], decision.values = T)
  )$decision.values
rocplot(-fitted, dat[-train, "y"], main = "Test Data")
fitted <- attributes(
    predict(svmfit.flex, dat[-train, ], decision.values = T)
  )$decision.values
rocplot(-fitted, dat[-train, "y"], add = T, col = "red")

Classificador de Vetores de Suporte com Múltiplas Classes

Multiclass Support Vector Classifier

Se nosso objetivo for utlizar uma SVM para classificação de dados em mais de \(2\) categorias, a função svm() realizará a classificação multiclasse usando a abordagem “Um-Contra-Um” (One-Versus-One). Exploramos essa configuração aqui gerando uma terceira classe de observações:

set.seed(1)
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))

Agora ajustamos uma SVM aos dados:

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

Agora examinamos o conjunto de dados Khan, da biblioteca ISLR2, que consiste em várias amostras de tecido correspondentes a quatro tipos distintos de tumores de células. Para cada amostra de tecido, medições de expressão gênica estão disponíveis. O conjunto de dados consiste em dados de treinamento, xtrain e ytrain, e dados de teste, xtest e ytest. Primeiro, examinamos a dimensão dos dados:

library(ISLR2)
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

Este conjunto de dados consiste em medições de expressão para \(2308\) genes. Os conjuntos de treinamento e teste consistem em \(63\) e \(20\) observações, respectivamente:

table(Khan$ytrain)

 1  2  3  4 
 8 23 12 20 
table(Khan$ytest)

1 2 3 4 
3 6 6 5 

Usaremos uma abordagem de vetores de suporte para prever o subtipo de câncer usando medições de expressão gênica. Neste conjunto de dados, há um número muito grande de características em relação ao número de observações. Isso sugere que devemos usar um kernel linear, porque a flexibilidade adicional do uso de um kernel não linear pode ser desnecessária…

dat <- data.frame(
    x = Khan$xtrain, 
    y = as.factor(Khan$ytrain)
  )
out <- svm(y ~ ., data = dat, kernel = "linear", 
    cost = 10)
summary(out)

Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  10 

Number of Support Vectors:  58

 ( 20 20 11 7 )


Number of Classes:  4 

Levels: 
 1 2 3 4
table(out$fitted, dat$y)
   
     1  2  3  4
  1  8  0  0  0
  2  0 23  0  0
  3  0  0 12  0
  4  0  0  0 20

Vemos que não há erros de treinamento. Na verdade, isso não é surpreendente, porque o grande número de variáveis em relação ao número de observações implica que é fácil encontrar hiperplanos que separam completamente as classes. Porém, como sempre, estamos mais interessados não no desempenho do classificador de vetores de suporte nas observações de treinamento, mas sim em seu desempenho nas observações de teste:

dat.te <- data.frame(
    x = Khan$xtest, 
    y = as.factor(Khan$ytest))
pred.te <- predict(out, newdata = dat.te)
table(pred.te, dat.te$y)
       
pred.te 1 2 3 4
      1 3 0 0 0
      2 0 6 2 0
      3 0 0 4 0
      4 0 0 0 5

Vemos que usar cost = 10 resulta em \(2\) erros no conjunto de teste nestes dados—não é um resultado ruim, mas poderiamos fazer melhor?

Regressão de Vetores de Suporte

Support Vector Regression

Por fim, vamos criar um modelo de regressão usando SVM. Lembre-se que a biblioteca e1071 também pode ser usada para realizar regressão de vetores de suporte. Para isso, basta que o vetor de resposta passado pela função svm() seja numérico em vez de um fator. No exemplo a seguir, vamos novamente trabalhar com o dados Boston, da biblioteca MASS, e prever o valor medv (ou seja, a mediana do valor das casas).

library(MASS)

Attaching package: 'MASS'
The following object is masked from 'package:ISLR2':

    Boston
data("Boston")
str(Boston)
'data.frame':   506 obs. of  14 variables:
 $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
 $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
 $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
 $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
 $ rm     : num  6.58 6.42 7.18 7 7.15 ...
 $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
 $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
 $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
 $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
 $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
 $ black  : num  397 397 393 395 397 ...
 $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
 $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
summary(Boston$medv)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   5.00   17.02   21.20   22.53   25.00   50.00 

É recomendável normalizar os dados antes:

# Normalização min-max:
normalize <- function(x) {
  return((x - min(x)) / (max(x) - min(x)))
}
Boston_norm <- as.data.frame(lapply(Boston, normalize))

Como sempre, primeiro, dividos os dados originais em um conjuto de treinamento e um de teste:

set.seed(123)
indice <- sample(1:nrow(Boston_norm), size = 0.8 * nrow(Boston_norm))
train_data <- Boston_norm[indice, ]
test_data <- Boston_norm[-indice, ]

Então, ajustamos nossa SVM com a função svm() (e um kernel radial, nesse caso). O argumento type = "eps-regression" especifica que estamos construindo um modelo de regressão com margem \(ε\) (epsilon):

modelo_svr <- svm(medv ~ ., data = train_data, type = "eps-regression", kernel = "radial")
summary(modelo_svr)

Call:
svm(formula = medv ~ ., data = train_data, type = "eps-regression", 
    kernel = "radial")


Parameters:
   SVM-Type:  eps-regression 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.07692308 
    epsilon:  0.1 


Number of Support Vectors:  276

Em seguida, examinamos o desempenho da nossa SVM:

predicoes <- predict(modelo_svr, test_data)
# Calcular métricas de desempenho
rmse <- sqrt(mean((predicoes - test_data$medv)^2))
cat("RMSE:", round(rmse, 4), "\n")
RMSE: 0.0781 
plot(test_data$medv, predicoes,
     xlab = "Valor Real", ylab = "Previsão",
     main = "SVR: Valor Real vs. Previsão", pch = 16, col = "blue")
abline(0, 1, col = "red", lwd = 2)  # Linha de 45°

Como vimos anteriormente, podemos utilizar a função tune() para encontrar os “melhores” valores dos hiperparâmetros da nossa SVM (nesse caso, de epsilon, cost e gamma) através de validação cruzada (note que gamma é um parâmetro da função kernel radial):

tuned <- tune(svm, medv ~ ., data = train_data,
              ranges = list(epsilon = c(0.01, 0.05, 0.10),
                            cost = c(1, 10, 100),
                            gamma = c(0.01, 0.05, 0.1)),
              type = "eps-regression", kernel = "radial")

summary(tuned)

Parameter tuning of 'svm':

- sampling method: 10-fold cross validation 

- best parameters:
 epsilon cost gamma
     0.1  100  0.05

- best performance: 0.006124718 

- Detailed performance results:
   epsilon cost gamma       error  dispersion
1     0.01    1  0.01 0.009723408 0.009791715
2     0.05    1  0.01 0.009819542 0.009880381
3     0.10    1  0.01 0.009733040 0.009736159
4     0.01   10  0.01 0.007069703 0.009039821
5     0.05   10  0.01 0.006935235 0.008850161
6     0.10   10  0.01 0.006910541 0.008641657
7     0.01  100  0.01 0.007489416 0.008906964
8     0.05  100  0.01 0.007331656 0.008666036
9     0.10  100  0.01 0.007035706 0.008401406
10    0.01    1  0.05 0.007730983 0.008961469
11    0.05    1  0.05 0.007629389 0.008876398
12    0.10    1  0.05 0.007559679 0.008890312
13    0.01   10  0.05 0.006530520 0.008127144
14    0.05   10  0.05 0.006363461 0.008020301
15    0.10   10  0.05 0.006484534 0.008158894
16    0.01  100  0.05 0.006781531 0.006157804
17    0.05  100  0.05 0.006479854 0.006246778
18    0.10  100  0.05 0.006124718 0.006323504
19    0.01    1  0.10 0.008377270 0.009422979
20    0.05    1  0.10 0.008311227 0.009393884
21    0.10    1  0.10 0.008233110 0.009527110
22    0.01   10  0.10 0.006892641 0.007685373
23    0.05   10  0.10 0.006751581 0.007438340
24    0.10   10  0.10 0.006645544 0.007269761
25    0.01  100  0.10 0.008634768 0.006213377
26    0.05  100  0.10 0.008118501 0.006073019
27    0.10  100  0.10 0.007603674 0.006211117
best_model <- tuned$best.model

# Avaliação no teste
pred_best <- predict(best_model, test_data)
rmse_best <- sqrt(mean((pred_best - test_data$medv)^2))
cat("RMSE:", round(rmse_best, 4), "\n")
RMSE: 0.0602 

Nosso novo modelo obteve um erro de teste mais baixo do que o modelo anterior, o que indica uma melhora—mas será que um modelo linear conseguiria um desempenho ainda melhor?

Exercício

Agora é sua vez de ajustar modelos SVM a um novo conjunto de dados. Você pode fazer isso utilizando seus próprios dados, dados disponíveis na internet ou em pacotes do R. Outra opção é utilizar os dados sobre desmatamento municipal em 2004, disponíveis em https://thaleswest.wixsite.com/home/tutorials (note que o site contém a descrição das variáveis deste conjunto de dados).

Siga os seguintes passos:

  • Primeiro, faça uma inspeção geral dos dados a serem utilizados.
  • Crie novas variáveis categóricas baseadas nos dados originais, se necessário.
  • Selecione uma variável dependente e ao menos \(15\) variáveis independentes para ajustar seu modelo SVM.
  • Realize alguns testes modificando os hiperparâmetros da função svm(), incluindo o tipo de função kernel.
  • Crie ao menos dois modelos SVM distintos e avalie suas performances. Qual modelo foi o melhor?