Laboratorio 9.6: Support Vector Machines
Paula Cazali
Fiabilidad
Se usará la librería e1071 para demostrar los vectores de clasificación y y SVM.
9.6.1 Support Vector Classifier
La función svm() puede ser usada para crear un modelo de clasificación con vectores de soporte cuando el parámetro kernel = 'lienar'. Esta función usa diferentes formulaciones de ciertas ecuaciones.
También se usa la función svm() para crear un support vector classifier para la función dada por el parámetro costo. Se comenzará generando las observaciones que perteneces a dos clases.
set.seed(1)
x <- matrix(rnorm(20*2), ncol = 2)
y <- c(rep(-1,10), rep(1,10))
x[y == 1] <- x[y == 1,] + 1Se comienza verificando cuales clases son linealmente separables
plot(x, col = (3-y))
Las clases no son separables. Por lo tanto se usara support vector classifier. Ahora se crea un data frame con la respuesta como factor.
dat <- data.frame(x = x, y = as.factor(y))
library(e1071)
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 10, scale = FALSE)El parámetro scale = FALSE indica a la función svm() que no escale cada variable para que tengan media igual a cero o desviación estandar igual a \(1\), dependiendo de la aplicación este parámetro puede tener TRUE o FALSE. Ahora se puede obtener el support vector classifier obtenido:
plot (svmfit, dat)
Notar que los dos argumentos graficar el SVM, es la salida de la función svm(). El límite entre las dos clases es lineal ya que se uso el argumento kernel = "linear".
svmfit$index[1] 1 2 5 7 14 16 17Se puede obtener información basica de la regresión del support vector classifier usando 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
gamma: 0.5
Number of Support Vectors: 7
( 4 3 )
Number of Classes: 2
Levels:
-1 1Esto nos dice que un kernel lineal fue usado con cost = 10 y que hay siete support vectors, tres en una clase y cuatro en la otra. Que pasa si usamos un valor más pequeño para el parámetro cost:
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 20Mientras se usa un valor más pequeño para ese parámetro, se obtiene un número mayor de support vectors, ya que el margen es más amplio. La función tune() crea \(10\) folds para realizar cross-validation en un set de modelos de interes.
set.seed(1)
tune.out <- tune(svm, y~ ., data = dat, kernel = "linear", ranges = list(cost = c(0.001, 0.01, 1, 5, 10, 100)))Ahora se puede acceder a los errores de cross-validation para cada uno de esos modelos usando summary():
summary(tune.out)
Parameter tuning of ‘svm’:
- sampling method: 10-fold cross validation
- best parameters:
cost
1
- best performance: 0.15
- Detailed performance results:
cost error dispersion
1 1e-03 0.70 0.4216370
2 1e-02 0.70 0.4216370
3 1e+00 0.15 0.2415229
4 5e+00 0.15 0.2415229
5 1e+01 0.15 0.2415229
6 1e+02 0.15 0.2415229Ahora se puede ver que al utilizar cost = 0.1 se obtiene el error más pequeño de cross-validation. Por lo tanto la función tune() puede obtener el mejor modelo, al cual se puede acceder:
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, 1, 5, 10,
100)), kernel = "linear")
Parameters:
SVM-Type: C-classification
SVM-Kernel: linear
cost: 1
gamma: 0.5
Number of Support Vectors: 11
( 6 5 )
Number of Classes: 2
Levels:
-1 1Usando la función predict() se puede predecir la etiqueta de la clase en un set de observaciones de prueba. Se generará un test:
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(predict = ypred, truth = testdat$y) truth
predict -1 1
-1 11 1
1 0 8Por lo que usando ese valor para el parámetro cost se pueden clasificar \(19\) observaciones de manera correcta. Qué hubiera pasado si se usara ahora \(cost = 0.01\)
svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 0.01, scale = FALSE)
ypred <- predict(svmfit, testdat)
table(predict = ypred, truth = testdat$y) truth
predict -1 1
-1 11 2
1 0 7En este caso una observación se pierde. Ahora considerando que dos clases son linealmente separables. Ahora se crearan dos clases que serán linealmente separables:
x[y == 1,] <- x[y == 1,] + 0.5
plot(x, col = (y + 5)/2, pch = 19)
Ahora las observaciones son casi separables linealmente.
dat <- data.frame(x = x, y = as.factor(y))
svmfit <- svm(y ~ ., data=dat, kernel ="linear", cost =1e+05)
summary(svmfit)
Call:
svm(formula = y ~ ., data = dat, kernel = "linear",
cost = 1e+05)
Parameters:
SVM-Type: C-classification
SVM-Kernel: linear
cost: 1e+05
gamma: 0.5
Number of Support Vectors: 3
( 1 2 )
Number of Classes: 2
Levels:
-1 1Ahora con un menor valor para 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
gamma: 0.5
Number of Support Vectors: 7
( 4 3 )
Number of Classes: 2
Levels:
-1 1plot(svmfit, dat)
9.6.2 Support Vector Machine
Para usar un kernel que no sea lineal también se usa la función svm() solo que ahora se usa un valor diferente para el parámetro kernel. Para un SVM polinomial usamos kernel = "polynomial" y para un kernel radial kernel = "radial", con el parámetro degree se especifica el grado para el kernel. Se genera data para aplicar un limite no lineal:
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))Graficamos:
plot(x, col = y)
Se separará la data en test y train:
train <- sample(200, 100)
svmfit <- svm(y ~ ., data = dat[train,], kernel = "radial", gamma = 1, cost = 1)
plot(svmfit, dat[train,])
La grafica muestra el limite no lineal creado.
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
gamma: 1
Number of Support Vectors: 37
( 17 20 )
Number of Classes: 2
Levels:
1 2Si se incrementa el valor para cost se pueden reducir el número de errores de training.
svmfit <- svm(y ~ ., data = dat[train,], kernel = "radial", gamma = 1, cost = 1e5)
plot(svmfit, dat[train,])
También se puede hacer cross-validation usando la función tune():
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 2
- best performance: 0.12
- Detailed performance results:
cost gamma error dispersion
1 1e-01 0.5 0.27 0.11595018
2 1e+00 0.5 0.13 0.08232726
3 1e+01 0.5 0.15 0.07071068
4 1e+02 0.5 0.17 0.08232726
5 1e+03 0.5 0.21 0.09944289
6 1e-01 1.0 0.25 0.13540064
7 1e+00 1.0 0.13 0.08232726
8 1e+01 1.0 0.16 0.06992059
9 1e+02 1.0 0.20 0.09428090
10 1e+03 1.0 0.20 0.08164966
11 1e-01 2.0 0.25 0.12692955
12 1e+00 2.0 0.12 0.09189366
13 1e+01 2.0 0.17 0.09486833
14 1e+02 2.0 0.19 0.09944289
15 1e+03 2.0 0.20 0.09428090
16 1e-01 3.0 0.27 0.11595018
17 1e+00 3.0 0.13 0.09486833
18 1e+01 3.0 0.18 0.10327956
19 1e+02 3.0 0.21 0.08755950
20 1e+03 3.0 0.22 0.10327956
21 1e-01 4.0 0.27 0.11595018
22 1e+00 4.0 0.15 0.10801234
23 1e+01 4.0 0.18 0.11352924
24 1e+02 4.0 0.21 0.08755950
25 1e+03 4.0 0.24 0.10749677Por lo tanto la mejor selección de parámetros se hace cuando \(cost = 1\) y \(gamma = 2\) Ahora se creará la matriz de confusión:
table(true = dat[-train,"y"], pred = predict(tune.out$best.model, newx = dat[-train,])) pred
true 1 2
1 56 21
2 18 59.6.3 ROC Curves
library(ROCR)Se creará una función para graficar la curva ROC dado un vector que contenga un puntaje numérico para cada observación:
rocplot <- function(pred, truth, ...){
predob <- prediction(pred, truth)
perf <- performance(predob, "tpr", "fpr")
plot(perf, ...)
}La función predict() devolverá los valores de salida del modelo:
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.valuesAhora se puede producir la gráfica:
par(mfrow = c(1,2))
rocplot(fitted, dat[train,"y"], main = " Training Data")
El SVM parece que produce predicciones exactas. Si aumentamos el valor de gamma se puede producir un modelo más flexible para mejorar la precision:
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"], col = "red")
Estas curvas ROC son sobre el training set. Ahora nos interesa en el nivel de exactitud de las predicciones del test set. Cuando se obtiene la curva ROC con \(gamma = 2\) parece que provee resultados más exactos.
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"], col = "red")
9.6.4 SVM with Multiple Classes
Si la respuesta es un factor que contiene más de dos niveles, entonces la función svm() puede realizar una clasificación multi-class usando el enfoque one-versus-one. Se genera una nueva clase de observaciones:
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))
Ahora se ajustará un SVM a la data:
svmfit <- svm(y ~ ., data = dat, kernel = "radial", cost = 10, gamma = 1)
plot(svmfit, dat)
La libreria e1071 puede ser usada para crear regresiones SV si el vector de respuesta que se para a la función svm() es numerica.
9.6.5 Application to Gene Expression Data
Ahora se examinara el dataset Khan que corresponden a tumores azules de una célula. El dataset consiste de data para train, xtrain y ytrain. y testing data xtest y ytest.
library(ISLR)
names(Khan)[1] "xtrain" "xtest" "ytrain" "ytest" dim(Khan$xtrain)[1] 63 2308dim(Khan$xtest)[1] 20 2308length(Khan$ytrain)[1] 63length(Khan$ytest)[1] 20El dataset consiste de mediciones para \(2,308\) genes. el training consiste de \(63\) observaciones y el testing consiste de \(20\).
table(Khan$ytrain)
1 2 3 4
8 23 12 20 table(Khan$ytest)
1 2 3 4
3 6 6 5 Se usará el enfoque de Support vector para predecir cancer usando la medición de los genes. En este dataset son muchas variables para pocas observaciones. Se usará un kernel lineal porque la flexibilidad adicional de usar un kernel polinomial o radial es innecesaria.
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
gamma: 0.0004332756
Number of Support Vectors: 58
( 20 20 11 7 )
Number of Classes: 4
Levels:
1 2 3 4table(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 20Por lo que se puede ver que no hay errores del training. Esto se puede deducir ya que el número de variables con las pocas observaciones pueden hacer que las clases se separen totalmente. Ahora se probará en el test set.
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 5Se puede ver que usando \(cost = 10\) se produce un error de \(2\) en estos datos.