SVM.Rmd

---
title: "SVM"
date: "2022/5/10"
output: pdf_document
---

```{r setup, include=FALSE}
library(tidyverse)
library(caret)
library(recipes)
library(vip)
library(pdp)
library(lime)
library(kernlab)

theme_set(theme_minimal() + theme(legend.position = "bottom"))

knitr::opts_chunk$set(echo = TRUE)
```


```{r}
hcv <- read.csv("HepatitisCdata.csv")[,-1]
#sum(is.na(hcv)) #31
hcv$Category <- case_when(hcv$Category == '0=Blood Donor' ~ 0,
                          hcv$Category == '0s=suspect Blood Donor' ~ 0,
                          TRUE ~ 1)
hcv$Sex <- factor(hcv$Sex, levels = c("f", "m"))


hcv.numeric <- hcv[,!names(hcv) %in% c("Category", "Sex")]
hcv.factor <- hcv[,!names(hcv) %in% colnames(hcv.numeric)]
hcv.factor <- hcv.factor[,!names(hcv.factor) %in% "Category"]

y <- hcv$Category
hcv$Category <- as.factor(ifelse(hcv$Category == 1, "patient", "donor"))
hcv$Sex <- factor(as.numeric(hcv$Sex)-1)
```


```{r}
set.seed(2022)
rowTrain <- createDataPartition(y = hcv$Category, p = 0.7, list = FALSE)
options(na.action='na.pass')

hcv.train.df <- hcv[rowTrain,]
hcv.test.df <- hcv[-rowTrain,]
```

```{r preprocess_with_impute}
preprocess_recipe = recipe(Category ~ ., data = hcv.train.df) %>%
    step_knnimpute(all_predictors(), neighbors = 5) %>%  
    step_BoxCox(all_numeric()) %>% 
    step_center(all_numeric()) %>% 
    step_scale(all_numeric())
```


```{r impute}
set.seed(2022)

prep <- prep(preprocess_recipe, training = hcv.train.df)
prep$template
# !! apply the preprocessing operations 
hcv.train.df <- bake(prep, new_data = hcv.train.df)
hcv.test.df <- bake(prep, new_data = hcv.test.df)

hcv.train.x <- model.matrix(~.-1,hcv[rowTrain,])
hcv.train.y <- hcv$Category[rowTrain]
hcv.test.x <- model.matrix(~.-1,hcv[-rowTrain,])
hcv.test.y <- hcv$Category[-rowTrain]
```


```{r}
ctrl <- trainControl(method = "repeatedcv",
                    repeats = 5)

ctrl_test = trainControl(method = "cv")
```


```{r}
library(parallel)
# Calculate the number of cores
num_cores <- detectCores() - 1
library(doParallel)
# create the cluster for caret to use
# CPU usage may go up to 100%
cl <- makePSOCKcluster(num_cores)
registerDoParallel(cl)

```


```{r}
set.seed(2022)
svml.fit = train(Category~.,
                 data = hcv.train.df,
                 method = "svmLinear",
                 metric = "Accuracy",
                 tuneGrid = data.frame(C = exp(seq(-2,2,len = 200))),
                 trControl = ctrl)

plot(svml.fit, highlight = TRUE, xTrans = log)
svml.fit$bestTune
#saveRDS(svml.fit, "svml.rds")

pred.svml = predict(svml.fit, newdata = hcv.test.df)
confusionMatrix(data = pred.svml,
                reference = hcv.test.df$Category)

```

```{r}
svmr.grid = expand.grid(C = exp(seq(-1,4,len = 50)),
                        sigma = exp(seq(-6,-2,len=50)))

set.seed(2022)
svmr.fit = train(Category~.,
                 data = hcv.train.df,
                 method = "svmRadialSigma",
                 tuneGrid = svmr.grid,
                 metric = "Accuracy",
                 trControl = ctrl)

myCol = rainbow(20)
myPar = list(superpose.symbol = list(col = myCol),
             superpose.line = list(col = myCol))

plot(svmr.fit, hightlight = TRUE, par.settings = myPar)

svmr.fit$bestTune
#saveRDS(svmr.fit, "svmr.rds")

pred.svmr = predict(svmr.fit, newdata = hcv.test.df)
confusionMatrix(data = pred.svmr,
                reference = hcv.test.df$Category)

```


```{r}
library(e1071)
test.fit = svm(Category~.,
               data = hcv.train.df,
               kernel = "radial",
               cost = svmr.fit$bestTune$C,
               gamma = svmr.fit$bestTune$sigma,
               #for pdp
               probability = TRUE)
predict(test.fit, newdata = hcv.test.df)

#classify
par(mfrow = c(1,2))
plot(test.fit, hcv.train.df,
     BIL ~ AST,
     slice = list(Age = 0, ALB = 0, ALP = 0, ALT = 0,
                  CHE = 0, CHOL = 0, CREA = 0, GGT = 1, PROT = 1),
     grid = 100)
plot(test.fit2, hcv.train.df,
     BIL ~ AST,
     slice = list(Age = 0, ALB = 0, ALP = 0, ALT = 0,
                  CHE = 0, CHOL = 0, CREA = 0, GGT = 1, PROT = 1),
     grid = 100)

#vip
w = t(test.fit$coefs) %*% test.fit$SV
w <- apply(w, 2, function(v){sqrt(sum(v^2))})
w <- sort(w)
print(w)

barplot(w, las = 2, horiz = TRUE, cex.names = .9, col = colorRampPalette(colors = c("cyan", "blue"))(13), xlim = c(0,25))

#pdp
features = c("AST", "BIL", "GGT", "PROT", "ALT", "ALP")

pdps <- lapply(features, function(x) {
  partial(test.fit, pred.var = x, which.class = 2,
          prob = TRUE, plot = TRUE, plot.engine = "ggplot2")
})
grid.arrange(grobs = pdps,  ncol = 3)
```


```{r}
library(rgl)
library(misc3d)

plot3d(hcv.train.df[c(5,6,7)], col = as.numeric(hcv.train.df$Category), size = 4)

# Get decision values for a new data grid
len = 10
newdat.list = lapply(hcv.train.df[c(5,6,7)], function(x) seq(min(x), max(x), len=len))
newdat      = expand.grid(newdat.list)

new = data.frame(Age = rep(0,len^3), Sex = c(rep(0,(len^3)/2), rep(1,(len^3)/2)), ALB = rep(0,len^3), ALP = rep(0, len^3), ALT = newdat$ALT, AST = newdat$AST, BIL = newdat$BIL, CHE = rep(0,len^3), CHOL = rep(0,len^3), CREA = rep(0,len^3), GGT = rep(1,len^3), PROT = rep(1.3, len^3))
new = as.tibble(new) %>% mutate(Sex = factor(Sex))

newdat.pred = predict(test.fit, newdata = new, decision.values = TRUE)
newdat.dv   = attr(newdat.pred, 'decision.values')
newdat.dv   = array(newdat.dv, dim=rep(len, 3))

# Fit/plot an isosurface to the decision boundary
contour3d(newdat.dv, level=0, x=newdat.list$ALT, y=newdat.list$AST, z=newdat.list$BIL, add=T, alpha = 0.5, color = "skyblue", smooth = 100)
```