Fiabilidad, Laboratorio 2

Laboratorio de Maquinas de Vectores de Soporte

Clasificador de vectores de soporte

La libreria e1071 contiene implementaciones para una serie de metodos estadisticos de aprendisaje. En particular, la función svm() se puede utilizar para realizar el ajuste:

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
x

            [,1]        [,2]
 [1,] -0.6264538  0.91897737
 [2,]  0.1836433  0.78213630
 [3,] -0.8356286  0.07456498
 [4,]  1.5952808 -1.98935170
 [5,]  0.3295078  0.61982575
 [6,] -0.8204684 -0.05612874
 [7,]  0.4874291 -0.15579551
 [8,]  0.7383247 -1.47075238
 [9,]  0.5757814 -0.47815006
[10,] -0.3053884  0.41794156
[11,]  2.5117812  2.35867955
[12,]  1.3898432  0.89721227
[13,]  0.3787594  1.38767161
[14,] -1.2146999  0.94619496
[15,]  2.1249309 -0.37705956
[16,]  0.9550664  0.58500544
[17,]  0.9838097  0.60571005
[18,]  1.9438362  0.94068660
[19,]  1.8212212  2.10002537
[20,]  1.5939013  1.76317575

Ahora revisemos si las clases son linealmente separables:

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

Ahora probaremos el ajuste del vector de clasificacion:

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

El argumento scale=FALSE le indica a svm() que no escale cada parametro para tener media 0 o desviacion standard 1, por lo tanto intentemos graficar esto:

plot(svmfit, dat)

Ahora determinemos los indices:

svmfit$index

[1]  1  2  5  7 14 16 17

Al ejecutar summary() podemos tener mejor informacion:

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

Esto nos confirma que tenemos un cost=10, que tenemos 7 vectores de soporte, 4 en una clase y 3 en la otra. Ahora probemos a utilizar un costo mas pequeno:

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

Al bajar el valor del costo vemosq ue obtenemos mas vectores de soporte. Ahora intenemos hacer un pequeno ajuste:

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)))
summary(tune.out)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:

- best performance: 0.1 

- Detailed performance results:

NA

El summary nos dice que cost=0.1 es el que menor error nos da.

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

Ahora utilicemos predict para validar nuestro modelo:

xtest <- matrix(rnorm(20*2), ncol=2)
ytest <- sample(c(-1,-1), 20, rep=TRUE)
xtest[ytest==1,] <- xtest[ytest==1,]+1
testdata <- data.frame(x=xtest, y=as.factor(ytest))
ypredict <- predict(bestmod, testdata)
table(predict=ypredict, truth=testdata$y)

       truth
predict -1
     -1 17
     1   3

Ahora utilicemos cost=0.01

svmfit <- svm(y~., data=dat, kernel="linear", cost=0.01,scale=FALSE)
ypredict <- predict(svmfit, testdata)
table(predict=ypredict, truth=testdata$y)

       truth
predict -1
     -1 18
     1   2

Ahora veamos la grafica de la separacion del hiperplano:

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

Vemos que las observaciones son a penas linealmente separables, ahora haremos el fit:

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

Ahora intentemoslo con un costo menor:

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)

Maquinas de Vectores de Soporte

Para hacer el ajuste de un modelo SVM en un kernel no linear, ahora veamos la siguiente prueba:

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))
plot(x, col=y)

Ahora hagamos una prueba de la funcion svm utilizando \(\gamma=1\):

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

Ahora vemos que que al no utilizar kernel linear obtenemos un limite radial, veamos el detalle:

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

Ahora veamos un incremente de costo:

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

Utilizando la funcion tune podemos hacer una validacion cruzada:

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

La mejor configucarion parece ser cost=1 y gamma=2, ahora hagamos unas predicciones:

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

    pred
true  1  2
   1 68  4
   2  3 25

Curvas ROC

La librerir ROCR nos permite graficar las areas bajo la curva, y para ello comenzaremos haciendo una pequena funcion:

library(ROCR)

Loading required package: gplots

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 observación. Sin embargo, también es posible obtener valores ajustados para cada observación, cuales son los puntajes numéricos utilizados para obtener las etiquetas de clase. Por ejemplo, en el caso de un clasificador de vectores de soporte, el valor ajustado para una observación \(X = (X_1, X_2, \dots, X_p)^T\) toma la forma \(\hat{\beta_0} + \hat{\beta_1}X_1 + \hat{\beta_2}X_2 +\dots + \hat{\beta_p}X_p\). Para una SVM con un kernel no lineal, la ecuación que produce el juste el valor se da en (9.23). En esencia, el signo del valor ajustado determina En qué lado del límite de decisión se encuentra la observación. por lo tanto, el relación entre el valor ajustado y la predicción de clase para un determinado la observación es simple: si el valor ajustado excede de cero, entonces la observación se asigna a una clase, y si es menor que cero que a la asignada otro. Para obtener los valores ajustados para un ajuste de modelo SVM dado, nosotros use decision.values decision.values=TRUE al ajustar svm(). Entonces la función de predicción() dará salida a los valores ajustados.

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

Ahora hagamos la grafica

SVM parece tener predicciones acertadas, haremos la prueba incrementando \(\gamma\), para producir resultados mas flexibles:

svmfit.flex <- svm(y~., data=dat[train,], kernel="radial", gamma=50, cost=1, decision.values =TRUE)
fitted <- attributes(predict(svmfit.flex, dat[train,], decision.values=TRUE))$decision.values
rocplot(fitted, dat[train,"y"], col="red")

Estas graficas estan generadas sobre la data de entrenamiento, ahora haremos la prueba con \(\gamma=2\) ya que este se ve que provee mejores resultados:

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

SVM con Multiples Clases

Si la respuesta es un factor que contiene más de dos niveles, la función svm() realizará una clasificación de múltiples clases utilizando el enfoque de uno contra uno. Exploramos ese entorno aquí generando 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) )

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

Aplicacion a Gene Expression Data

Ahora analizaremos el dataset Khan

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

table(Khan$ytrain)

 1  2  3  4 
 8 23 12 20 

table(Khan$ytest)

1 2 3 4 
3 6 6 5 

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 relación con la cantidad de observaciones implica que es fácil encontrar hiperplanos que separen completamente las clases. Lo que más nos interesa no es el desempeño del clasificador de vectores de soporte en las observaciones de entrenamiento, sino 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

Ahora vemos que cost=10 arroja dos errores en la data.