Laboratorio 9.6 - Maquinas de vectores de soporte

Cesar Tinoco - 13003387

Se utiliza la libreria e1071 en R para demostrar el clasificador de vectores de soporte y el SVM. Otra opción es la libreria LiblineaR, que es útil para problemas lineales muy grandes.

9.6.1 Clasificador de vectores de soporte

La libreria e1071 contiene implementaciones para varios metodos de aprendizaje estadistico, en particular la función svm() un clasificador de vectores de soporte cuando se utilice el argumento kernel=“linear”. Un argumento de costos nos permite especificar el costo de violación al margen, cuando el argumento de costo es pequeño, entonces los margenes serán anchos. Ahora usamos la funcion svm () para ajustar el clasificador de vectores de soporte para un Valor dado del parametro coste.Comenzamos generando las observaciones, que pertenecen 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,] + 1

Se verifica si las clases son linealmente separables.

plot(x, col =(3-y))

Se pudo comprobar que no son, luego se ajusta el clasificador de vectores de soporte. Tenga en cuenta que para que la funcion svm() realice una clasificacion se debe codificar la respuesta como una variable de factor. A continuacion se crea un dataframe con la respuesta codificada como un factor.

library (e1071)

package 㤼㸱e1071㤼㸲 was built under R version 3.5.3

dat=data.frame(x=x, y=as.factor (y))
svmfit =svm(y∼., data=dat , kernel ="linear", cost =10,scale =FALSE )

El argumento scale = FALSE le dice a la funcion svm() que no escale cada caracteristica para que tenga una media de cero o una desviacion estandar de uno; Dependiendo de la aplicacion, uno podria preferir usar scale = TRUE

dat=data.frame(x=x, y=as.factor (y))
svmfit =svm(y∼., data=dat , kernel ="linear", cost =10,scale = TRUE )

Ahora podemos trazar el clasificador de vectores de soporte obtenido.

plot(svmfit , dat)

Hay que tener en cuenta que los argumentos de la funcion plot.svm() son la salida de la llamada a svm().La region del espacio de la caracteristica que se asignara a la clase −1 se muestra en azul claro, y la region que se asignara a la clase +1 se muestra en purpura. El límite de decision entre las dos clases es lineal (porque usamos el argumento kernel = “linear”), aunque debido a la forma en que se implementa la funcion de trazado en esta libreria, el limite de decision parece algo irregular en el trazado. Vemos que en este caso solo una observacion está mal clasificada.Los vectores de soporte se trazan como cruces y Las observaciones restantes se trazan como círculos; Vemos aquí que hay siete vectores de apoyo.

Podemos determinar las identidades de los vectores de la siguiente manera:

svmfit$index

[1]  1  2  5  7 14 16 17

Podemos obtener información sobre el ajuste del clasificador de vectores de soporte usando summary()

summary (svmfit)

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


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

Number of Support Vectors:  7

 ( 4 3 )


Number of Classes:  2 

Levels: 
 -1 1

De lo anterior podemos decir que se uso un kernel lineal con el costo = 10, y que habia siete vectores de soporte, cuatro en una clase y tres en la otra.

¿Qué pasaría si en cambio usáramos un valor más pequeño del parámetro de costo?

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

Como se esta utilizando un valor mas pequeño del parametro de costo se obtiene un numero mayor de vectores de soporte debido a que el ancho del margen es mas ancho. svm() no genera explicitamente los coeficientes del limite de decision lineal que se obtiene cuando el clasificador de vectores esta ajustado. La funcion tune() permite realizar cruces de validacion, en forma predeterminada realiza una validacion cruzada de diez veces en un conjunto de modelos de interes. Para utilizar esta funcion se transfiere informacion relevante del conjunto de modelos que se esta considerando. A continuacion realizaremos una comparacion entre SVM y un Kernel Lineal utilizando un rango de valores del parametro de costo.

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 acceder facilmente a los errores de validacion cruzada para cada uno de estos modelos usando summary():

summary (tune.out)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:

- best performance: 0.1 

- Detailed performance results:

NA

Vemos que el costo = 0.1 es la tasa de error de validacion cruzada mas baja. La funcion tune() almacena el mejor modelo obtenido, al cual se puede acceder de la siguiente manera:

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

La función predict() se puede usar para predecir la etiqueta de clase en un conjunto de observaciones de prueba, en cualquier valor dado del parametro de costo. Ahora generaremos un conjunto de datos de prueba.

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

Ahora predecimos las etiquetas de clase de estas observaciones de prueba. Aquí utilizamos el mejor modelo obtenido mediante validación cruzada para hacer predicciones.

ypred=predict (bestmod ,testdat )
table(predict =ypred , truth= testdat$y )

       truth
predict -1  1
     -1 11  1
     1   0  8

Por lo tanto, con este valor de costo, hay 19 observaciones de prueba que estan correctamente clasificadas. ¿Qué pasaría si hubiésemos utilizado en su lugar el costo = 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  2
     1   0  7

En este caso, una observación adicional está mal clasificada. Ahora considere una situación en la que las dos clases son linealmente separables. Luego podemos encontrar un hiperplano de separación usando la función svm(). Primero separamos las dos clases en nuestros datos simulados para que sean linealmente separables:

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

Ahora las observaciones son apenas linealmente separables. Ajustamos el clasificador de vectores de soporte y trazamos el hiperplano resultante, utilizando un valor de costo muy grande para que no se clasifiquen erróneamente las observaciones

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

No se cometieron errores de entrenamiento y solo se usaron tres vectores de soporte. Sin embargo, podemos ver en la figura que el margen es muy estrecho (porque las observaciones que no son vectores de soporte, indicadas como círculos, están muy cerca del límite de decision). Parece probable que este modelo se desempeñe mal en los datos de prueba. Ahora intentamos un menor valor de costo:

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 )

Usando cost = 1, clasificamos erroneamente una observacion de entrenamiento, pero tambien obtenemos un margen mucho mas amplio y utilizamos siete vectores de soporte. Parece probable que este modelo tenga un mejor desempeño en los datos de prueba que el modelo con cost = 1e5.

9.6.2 Maquinas de vectores soporte

Para ajustar un SVM usando un kernel no lineal, una vez mas usamos la funcion svm(). Sin embargo, ahora usamos un valor diferente del parámetro kernel. Para ajustar una SVM con un núcleo polinomial usamos kernel = “polinomial”, y para ajustar una SVM con un núcleo radial usamos kernel = “radial”. En el primer caso, también usamos el argumento de grado para especificar un grado para el núcleo polinomial (esto es d en (9.22)), y en el último caso usamos gamma para especificar un valor de γ para el núcleo de base radial (9.24) . Primero generamos algunos datos con un límite de clase 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))

El trazado de los datos deja claro que el limite de la clase es no lineal:

plot(x, col=y)

Los datos se dividen aleatoriamente en grupos de train y test. Entonces encajamos los datos del train utilizando la función svm () con un núcleo radial y γ = 1:

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

La grafica muestra que el SVM resultante tiene un limite no lineal. La funcion de summary() se puede utilizar para obtener alguna informacion sobre el ajuste 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:  35

 ( 15 20 )


Number of Classes:  2 

Levels: 
 1 2

Podemos ver en la figura que hay un buen numero de errores de entrenamiento en este ajuste de SVM. Si aumentamos el valor del costo, podemos reducir el número de errores de capacitacion. Sin embargo, esto viene al precio de un limite de decision mas irregular que parece estar en riesgo de sobreajustar los datos.

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

Podemos realizar una validación cruzada utilizando tune() para seleccionar la mejor opción de γ y el costo de una SVM con un núcleo 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:

- best performance: 0.09 

- Detailed performance results:

NA

Por lo tanto, la mejor eleccion de parametros implica cost = 1 y gamma = 2. Podemos ver las predicciones del conjunto de pruebas para este modelo aplicando la funcion predict() a los datos. Tener en cuenta que para hacer esto, subcontratamos los datos del marco de datos utilizando -train como un conjunto de indices.

table(true=dat[-train ,"y"], pred=predict (tune.out$best.model ,
newx=dat[-train ,]))

    pred
true  1  2
   1 55 23
   2 16  6

39% de las observaciones de prueba estan mal clasificadas por este SVM.

9.6.3 Curvas ROC

El paquete ROCR se puede utilizar para producir curvas ROC. Primero escribimos una función corta para trazar una curva ROC dado un vector que contiene una puntuacion numerica para cada observacion, pred y un vector que contiene la etiqueta de la clase para cada observación, truth.

library (ROCR)

package 㤼㸱ROCR㤼㸲 was built under R version 3.5.3Loading required package: gplots
package 㤼㸱gplots㤼㸲 was built under R version 3.5.3
Attaching package: 㤼㸱gplots㤼㸲

The following object is masked from 㤼㸱package:stats㤼㸲:

    lowess

rocplot =function (pred , truth , ...){
   predob = prediction (pred , truth )
   perf = performance (predob , "tpr", "fpr")
   plot(perf ,...)}

Los SVM y los clasificadores de vectores dan salida a las etiquetas de clase para cada observacion y tambien es posible obtener valores ajustados para cada observacion, que son los puntajes numericos utilizados para obtener las etiquetas de clase. Para una SVM con un kernel no lineal. En esencia, el signo del valor ajustado determina en que lado del limite de decision se encuentra la observacion. Por lo tanto, la relacion entre el valor ajustado y la prediccion de clase para una observacion dada es simple: si el valor ajustado excede de cero, la observacion se asigna a una clase y si es menor que cero se asigna a la otra. Para obtener los valores ajustados para un ajuste de modelo SVM dado, usamos decision.values = TRUE cuando ajustamos svm(). Entonces la función de predict() dará salida a los valores ajustados.

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

Ahora podemos producir ROC plot.

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

SVM parece estar produciendo predicciones precisas. Al aumentar γ podemos producir un ajuste mas flexible y generar mas mejoras en 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")

Sin embargo, estas curvas ROC estan todas en los datos de entrenamiento. Interesa mas el nivel de precision de prediccion en los datos de prueba. Cuando calculamos las curvas ROC en los datos de prueba, el modelo con γ = 2 parece proporcionar los resultados más 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 ")

9.6.4 SVM con multiples clases

Si la respuesta es un factor que contiene mas de dos niveles, la función svm() realizara una clasificacion de multiples clases utilizando el enfoque de uno contra uno. Se genera una tercera 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 ajustamos un SVM a los datos:

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

La libreria e1071 tambien se puede usar para realizar una regresion de vectores de soporte, si el vector de respuesta que se pasa a svm() es numerico en lugar de un factor.

9.6.5 Aplicacion a los datos de expresion genica

Ahora examinamos el conjunto de datos de Khan, que consiste en una serie de muestras de tejido que corresponden a cuatro tipos distintos de tumores pequeños redondos de células azules. Para cada muestra de tejido, las mediciones de expresion genica estan disponibles. El conjunto de datos consta de datos de entrenamiento, xtrain y ytrain, y datos de prueba, xtest y ytest. Examinamos la dimension de los datos:

library (ISLR)

package 㤼㸱ISLR㤼㸲 was built under R version 3.5.3

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 datos consiste en mediciones de expresion para 2,308 genes. Los conjuntos de entrenamiento y prueba constan de 63 y 20 observaciones de forma espectativa.

table(Khan$ytrain )

 1  2  3  4 
 8 23 12 20 

table(Khan$ytest )

1 2 3 4 
3 6 6 5 

Utilizaremos un enfoque de vectores de apoyo para predecir el subtipo de cancer mediante mediciones de expresion genica. En este conjunto de datos, hay una gran cantidad de caracteristicas en relacion con la cantidad de observaciones. Esto sugiere que deberiamos usar un nucleo lineal, porque la flexibilidad adicional que resultara del uso de un nucleo 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 

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 no hay errores de entrenamiento. De hecho, esto no es sorprendente, debido a que la gran cantidad de variables en relacion con la cantidad de observaciones implica que es facil encontrar hiperplanos que separan las clases por completo. Estamos mas interesados no en el desempeño del clasificador de vectores de soporte en las observaciones de entrenamiento, sino en su desempeño en las observaciones de prueba.

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 produce dos errores de conjunto de prueba en estos datos.