5.4 SVM en R

Hay varios paquetes que implementan este procedimiento (p. ej. e1071, Meyer et al., 2020; svmpath, ver Hastie et al., 2004), aunque se considera que el más completo es kernlab (Karatzoglou et al., 2004).

La función principal del paquete kernlab es ksvm() y se suelen utilizar los siguientes argumentos:

ksvm(formula, data, scaled = TRUE, type, kernel ="rbfdot", kpar = "automatic",
     C = 1, epsilon = 0.1, prob.model = FALSE, class.weights, cross = 0)

formula y data (opcional): permiten especificar la respuesta y las variables predictoras de la forma habitual (p. ej. respuesta ~ .; también admite matrices).
scaled: vector lógico indicando qué predictores serán reescalados; por defecto, se reescalan todas las variables no binarias (y se almacenan los valores empleados para ser usados en posteriores predicciones).
type (opcional): cadena de texto que permite seleccionar los distintos métodos de clasificación, de regresión o de detección de atípicos implementados (ver ?ksvm); por defecto, se establece a partir del tipo de la respuesta: "C-svc", clasificación con parámetro de coste, si es un factor, y "eps-svr", regresión épsilon, si la respuesta es numérica.
kernel: función núcleo. Puede ser una función definida por el usuario o una cadena de texto que especifique una de las implementadas en el paquete (ver ?kernels); por defecto, "rbfdot", kernel radial gaussiano.
kpar: lista con los hiperparámetros del núcleo. En el caso de "rbfdot", además de una lista con un único componente "sigma" (inversa de la ventana), puede ser "automatic" (valor por defecto) e internamente emplea la función sigest() para seleccionar un valor “adecuado”.
C: hiperparámetro \(C\) que especifica el coste de la violación de las restricciones; 1 por defecto.
epsilon: hiperparámetro \(\epsilon\) empleado en la función de pérdidas de los métodos de regresión; 0.1 por defecto.
prob.model: si se establece a TRUE (por defecto es FALSE), se emplean los resultados de la clasificación para ajustar un modelo para estimar las probabilidades (y se podrán calcular con el método predict()).
class.weights: vector (con las clases como nombres) con los pesos asociados a las observaciones de cada clase (por defecto, 1). Se puede entender como el coste de una mala clasificación en cada clase y podría ser de utilidad también para clases desbalanceadas.

cross: número de grupos para validación cruzada; 0 por defecto (no se hace validación cruzada). Si se asigna un valor mayor que 1, se realizará validación cruzada y se devolverá el error en la componente @cross (se puede acceder con la función cross(); y se puede emplear para seleccionar hiperparámetros).

Como ejemplo consideraremos el problema de clasificación con los datos de calidad de vino:

library(mpae)
data("winetaste")
# Partición de los datos
set.seed(1)
df <- winetaste
nobs <- nrow(df)
itrain <- sample(nobs, 0.8 * nobs)
train <- df[itrain, ]
test <- df[-itrain, ]
# Entrenamiento
library(kernlab)
set.seed(1) 
svm <- ksvm(taste ~ ., data = train, kernel = "rbfdot", prob.model = TRUE)
svm

## Support Vector Machine object of class "ksvm" 
## 
## SV type: C-svc  (classification) 
##  parameter : cost C = 1 
## 
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.0751133799772488 
## 
## Number of Support Vectors : 594 
## 
## Objective Function Value : -494.14 
## Training error : 0.198 
## Probability model included.

Podemos evaluar la precisión en la muestra de test empleando el procedimiento habitual:

pred <- predict(svm, newdata = test)
caret::confusionMatrix(pred, test$taste)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction good bad
##       good  147  45
##       bad    19  39
##                                         
##                Accuracy : 0.744         
##                  95% CI : (0.685, 0.797)
##     No Information Rate : 0.664         
##     P-Value [Acc > NIR] : 0.00389       
##                                         
##                   Kappa : 0.379         
##                                         
##  Mcnemar's Test P-Value : 0.00178       
##                                         
##             Sensitivity : 0.886         
##             Specificity : 0.464         
##          Pos Pred Value : 0.766         
##          Neg Pred Value : 0.672         
##              Prevalence : 0.664         
##          Detection Rate : 0.588         
##    Detection Prevalence : 0.768         
##       Balanced Accuracy : 0.675         
##                                         
##        'Positive' Class : good          
##

Para obtener las estimaciones de las probabilidades, habría que establecer type = "probabilities" al predecir (devolverá una matriz con columnas correspondientes a los niveles)⁵⁰:

p.est <- predict(svm, newdata = test, type = "probabilities")
head(p.est)

##         good     bad
## [1,] 0.47619 0.52381
## [2,] 0.70893 0.29107
## [3,] 0.88935 0.11065
## [4,] 0.84240 0.15760
## [5,] 0.66409 0.33591
## [6,] 0.36055 0.63945

Ejercicio 5.1 En el ejemplo anterior, el error de clasificación en la categoría bad es mucho mayor que en la otra categoría (la especificidad es 0.4643). Esto podría ser debido a que las clases están desbalanceadas y el modelo trata de clasificar mejor la clase mayoritaria. Podríamos tratar de mejorar la especificidad empleando el argumento class.weights de la función ksvm(). Por ejemplo, de forma que el coste de una mala clasificación en la categoría bad sea el doble que el coste en la categoría good⁵¹. De esta forma sería de esperar que se clasifique mejor la clase minoritaria, bad (que aumente la especificidad), a expensas de una disminución en la sensibilidad (para la clase mayoritaria good). Se esperaría también una mejora en la precisión balanceada, aunque con una reducción en la precisión. Repite el ejemplo anterior empleando el argumento class.weights para mejorar la especificidad.

Las SVM también están implementadas en caret, en múltiples métodos. Uno de los más empleados es "svmRadial" (equivalente a la clasificación anterior con núcleos radiales gaussianos) y considera como hiperparámetros:

library(caret)
# names(getModelInfo("svm")) # 17 métodos
modelLookup("svmRadial")

##       model parameter label forReg forClass probModel
## 1 svmRadial     sigma Sigma   TRUE     TRUE      TRUE
## 2 svmRadial         C  Cost   TRUE     TRUE      TRUE

En este caso, la función train() por defecto evaluará únicamente tres valores del hiperparámetro C = c(0.25, 0.5, 1) y fijará el valor de sigma. Alternativamente, podríamos establecer la rejilla de búsqueda. Por ejemplo, fijamos sigma al valor por defecto⁵² de la función ksvm() e incrementamos los valores del hiperparámetro de coste:

tuneGrid <- data.frame(sigma = kernelf(svm)@kpar$sigma, # Emplea clases S4
                       C = c(0.5, 1, 5))
set.seed(1)
caret.svm <- train(taste ~ ., data = train, method = "svmRadial", 
                   preProcess = c("center", "scale"),
                   trControl = trainControl(method = "cv", number = 5),
                   tuneGrid = tuneGrid, prob.model = TRUE)
caret.svm

## Support Vector Machines with Radial Basis Function Kernel 
## 
## 1000 samples
##   11 predictor
##    2 classes: 'good', 'bad' 
## 
## Pre-processing: centered (11), scaled (11) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 800, 801, 800, 800, 799 
## Resampling results across tuning parameters:
## 
##   C    Accuracy  Kappa  
##   0.5  0.75495   0.42052
##   1.0  0.75993   0.42975
##   5.0  0.75494   0.41922
## 
## Tuning parameter 'sigma' was held constant at a value of 0.075113
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.075113 and C = 1.

En este caso, el modelo obtenido es el mismo que en el ejemplo anterior (se seleccionaron los mismos valores de los hiperparámetros).

Ejercicio 5.2 Repite el ajuste anterior realizando una búsqueda de ambos hiperparámetros para tratar de mejorar la clasificación. Para el hiperparámetro sigma se podría considerar como referencia el valor seleccionado automáticamente por la función ksvm(). Por ejemplo, incluyendo además la mitad y el doble de ese valor: sigma = with(kernelf(svm)@kpar, c(sigma/2, sigma, 2*sigma).

Ejercicio 5.3 Para seleccionar los hiperparámetros en un problema de clasificación, caret utiliza como criterio por defecto la precisión de las predicciones. En la Sección 1.3.5 se mostraron criterios alternativos que podrían resultar de interés en ciertos casos. Por ejemplo, para emplear el área bajo la curva ROC (AUC), en primer lugar necesitaríamos añadir classProbs = TRUE en trainControl(), ya que esta medida precisa de las estimaciones de las probabilidades de cada clase, que no se calculan por defecto. En segundo lugar, habría que cambiar la función que calcula los distintos criterios de optimalidad en la llamada a trainControl(). Estableciendo summaryFunction = twoClassSummary se calcularían medidas específicas para problemas de dos clases: el área bajo la curva ROC, la sensibilidad y la especificidad (en lugar de la precisión y el valor de Kappa). Finalmente, habría que incluir metric = "ROC" en la llamada a train() para establecer el AUC como criterio de selección de hiperparámetros.

Repite el ejemplo anterior seleccionando los hiperparámetros de forma que se maximice el área bajo la curva ROC.

Ejercicio 5.4 Continuando con el conjunto de datos mpae::bfan empleado en ejercicios de capítulos anteriores, particiona los datos y clasifica los individuos según su nivel de grasa corporal (bfan):

Empleando la función ksvm() del paquete kernlab con las opciones por defecto.
Empleando el método "svmRadial" del paquete caret, seleccionando los valores óptimos de los hiperparámetros mediante validación cruzada con 5 grupos y considerando las posibles combinaciones de C = c(0.5, 1, 5) y sigma = c(0.5, 1, 2)*sigma0, siendo sigma0 el valor de este parámetro seleccionado en el apartado anterior.
Evalúa la precisión de las predicciones de ambos modelos en la muestra de test y compara los resultados.

Bibliografía

Hastie, T., Rosset, S., Tibshirani, R., y Zhu, J. (2004). The entire regularization path for the support vector machine. Journal of Machine Learning Research, 5, 1391-1415.

Karatzoglou, A., Smola, A., Hornik, K., y Zeileis, A. (2004). kernlab - An S4 Package for Kernel Methods in R. Journal of Statistical Software, 11(9), 1-20. https://doi.org/10.18637/jss.v011.i09

Meyer, D., Dimitriadou, E., Hornik, K., Weingessel, A., y Leisch, F. (2020). e1071: Misc Functions of the Department of Statistics, Probability Theory Group (Formerly: E1071), TU Wien. https://cran.r-project.org/package=e1071

Otras opciones son "votes" y "decision" para obtener matrices con el número de votos o los valores de \(m(\mathbf{x})\).↩︎
Una alternativa similar, que se suele emplear cuando las clases están desbalanceadas, es ponderar las observaciones por un peso inversamente proporcional a la frecuencia de cada clase.↩︎
La función ksvm(), por defecto, selecciona sigma = mean(sigest(taste ~ ., data = train)[-2]), aunque hay que tener en cuenta que el resultado de la función sigest() depende de la semilla.↩︎