#####################################################################################################
##
##                                                SVM
##
## E. Flachaire, 13/12/2021
##
#####################################################################################################


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))
plot(x[,2], x[,1], pch=16, col=y*2)

# SVM with 2nd degree polynomial kernel and fixed cost of contraints violation
library(e1071)
dat <- data.frame(x=x, y=as.factor(y))
svmfit <- svm(y~., data=dat, kernel="polynomial", cost=1, degree=2)
plot(svmfit, dat, grid=200)

# Select the cost parameter by 10-fold CV among several values
tune.out <- tune(svm, y~., data=dat, kernel="polynomial", degree=2,
                 ranges=list(cost=c(.1,1,10,100)))
plot(tune.out$best.model, dat, grid=200)
summary(tune.out)

# SVM with radial kernel
tune.out <- tune(svm, y~., data=dat, kernel="radial", 
                 ranges=list(cost=c(.1,1,10,100),
                             gamma=c(.5,1,2,3,4)))
plot(tune.out$best.model, dat, grid=200)

# ROC curves

library(ROCR)

#' Short function to plot a ROC curve
#' @param pred Vector or predicted values
#' @param truth Vector of true (observed) values
#' @param ... Other arguments passed on to methods
rocplot <- function(pred, truth, ...){
  predob <- prediction(pred, truth)
  perf <- performance(predob, "tpr", "fpr")
  plot(perf, ...)
}


# Fit a SVM with radial kernel and plot a ROC curve
set.seed(1)
train <- sample(200,100)
train <- sort(train, decreasing=TRUE) # to avoid reverse ROC
svmfit <- svm(y~., data=dat[train,], kernel = "radial", cost=1, gamma=0.5)
fit <- attributes(predict(svmfit, dat[-train,], decision.values=TRUE))$decision.values
rocplot(fit, dat[-train, "y"], main = "Test Data", col="red")

# Logit model for comparison
lgt <- glm(y~., data=dat[train,], family=binomial(link="logit"))
fit <- predict(lgt, dat[-train,], type="response")
par(new=TRUE)
rocplot(fit, dat[-train, "y"], col="green")