rm(list = ls()); ls()
## character(0)
library(e1071)
library(ROCR)
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## 载入程辑包：'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(ROSE)
## Loaded ROSE 0.0-4
library(DMwR)
## 载入需要的程辑包：lattice
## 载入需要的程辑包：grid
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
###############################input data_1
dir_path <- "C:\\Users\\liyix\\OneDrive\\Desktop\\data_for_model\\"
dir_path_name <- dir(dir_path,pattern = ".*.csv",full.names = T)
dir_path_name
## [1] "C:\\Users\\liyix\\OneDrive\\Desktop\\data_for_model\\2021-12-10-nb_upsample_probability.csv"
## [2] "C:\\Users\\liyix\\OneDrive\\Desktop\\data_for_model\\2021-12-10-nnet_upsample_probability.csv"
## [3] "C:\\Users\\liyix\\OneDrive\\Desktop\\data_for_model\\test_data.csv"
## [4] "C:\\Users\\liyix\\OneDrive\\Desktop\\data_for_model\\training_data.csv"
###############################merge data
data_train <- read.csv(grep("training_data.csv",dir_path_name,value = T),header = T,stringsAsFactors = F)
dim(data_train) #[1] 2357 148
## [1] 2357 148
#View(data_train)
data_test <- read.csv(grep("test_data.csv",dir_path_name,value = T),header = T,stringsAsFactors = F)
#dim(data_test) #[1] 1009 150
#table(data_train$Freq) # 0 1 2057 300
#table(data_test$Freq) # 0 1 894 115
#View(head(data_test))
#setdiff(colnames(data_test),colnames(data_train)) #[1] "Drug" "Virus"
#setdiff(colnames(data_train),colnames(data_test)) #character(0)
#colnames(data_train)[ncol(data_train)] #[1] "Freq"
colnames(data_train)[ncol(data_train)] <- "num"
data_train <- data.frame(sapply(data_train, as.numeric))
data_train$num <- as.factor(data_train$num)
#str(data_train$num)
############################################################################test_cross_validation
#dim(data_train) #[1] 2356 148
#View(head(data_train))
#colnames(data_train)
#table(data_train$num) # 0 1 2057 300
testdata_ori <- data_train[sample(nrow(data_train),size = floor(nrow(data_train)*0.3),replace = F),]
traindata_ori <- data_train[-match(row.names(testdata_ori),row.names(data_train)),]
#intersect(row.names(testdata_ori),row.names(traindata_ori)) #character(0)
#table(testdata_ori$num) # 0 1 625 81
#table(traindata_ori$num) # 0 1 1444 206
#dim(testdata_ori);dim(traindata_ori) #[1] 706 148 [1] 1650 148
##model
#newdata_tune_svm <- tune.svm(num ~ .,data = train_data, cost = 10^(-3:3),
# gamma = c(0.01,0.05,0.1,0.5,1), tunecontrol =
# tune.control(cross = 5), scale = F)
#newdata_tune_svm$best.parameters$gamma <- 0.005
#newdata_tune_svm$best.parameters$cost <- 10
#################################################probability
traindata_ori <- ovun.sample(num ~ ., data = traindata_ori, method = "under", seed = 1)$data
fit_nb_ori <- naiveBayes(num ~ .,data = traindata_ori,laplace = 1)
#########################################################
#### auc_rf
n_col_newdata <- ncol(data_train)
pre_auc_nb_ori <- predict(fit_nb_ori, testdata_ori[,-n_col_newdata], type = 'raw')
pred_auc_nb_ori <- prediction(pre_auc_nb_ori[,2],testdata_ori[,n_col_newdata])
auc_nb_ori <- performance(pred_auc_nb_ori,"auc")@y.values[[1]]
auc_nb_ori #[1] 0.9793877
## [1] 0.7220337
pre_probability_ori <- data.frame(pre_auc_nb_ori)
pre_classification_ori <- predict(fit_nb_ori, testdata_ori[,-ncol(data_test)], type="class")
#class(pre_classification_ori)
pre_probability_ori$num_1 <- as.character(pre_classification_ori)
pre_probability_ori$num <- testdata_ori$num
#head(pre_probability_ori)
#table(pre_probability_ori$num,pre_probability_ori$num_1)
modelroc_nb_ori <- pROC::roc(testdata_ori[,n_col_newdata],as.numeric(pre_auc_nb_ori[,2]))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
cutoff_nb_ori <- coords(modelroc_nb_ori, "best", ret = "threshold") #[1] 0.265
predict.results_nb_ori <- ifelse(pre_auc_nb_ori[,2] > as.numeric(cutoff_nb_ori[1]),"1","0")
freq_default_nb_ori <- table(predict.results_nb_ori,testdata_ori[,n_col_newdata])
Sensitivity_nb_ori <- freq_default_nb_ori[1,1]/sum(freq_default_nb_ori[,1])
Specificity_nb_ori <- freq_default_nb_ori[2,2]/sum(freq_default_nb_ori[,2])
accuracy_nb_ori <- mean(testdata_ori[,n_col_newdata] == predict.results_nb_ori)
Balanced_accuracy_nb_ori <- (Sensitivity_nb_ori+Specificity_nb_ori)/2
Balanced_accuracy_nb_ori #[1] 0.954963
## [1] 0.7145903
library(caret)
## 载入需要的程辑包：ggplot2
#confusionMatrix(factor(predict.results_nb_ori),factor(testdata_ori[,n_col_newdata]))
#confusionMatrix(factor(pre_probability_ori$num_1),factor(testdata_ori[,n_col_newdata]))
#####0 1
#0 1350 15
#1 116 13
#################################################################################################
##nb
data_train <- ovun.sample(num ~ ., data = data_train, "under", seed = 1)$data
#dim(data_train) #[1] 4129 148
#################################################probability
#fit_probability <- svm(num ~ .,data = data_train, scale = F, probability=TRUE )
#dim(data_train) #[1] 4981 116
#table(data_train$num) #4515 466
fit_probability <- naiveBayes(num ~ .,data = data_train, laplace = 1)
#### AUC
#names(data_test)[ncol(data_test)]
#dim(data_test) #[1] 1009 150
#colnames(data_test)
pre_auc_probability <- predict(fit_probability, data_test[,-c(ncol(data_test),ncol(data_test)-1, ncol(data_test)-2)], type = 'raw')
pre_probability <- data.frame(pre_auc_probability)
#head(pre_probability)
########################################
colnames(pre_probability) <- gsub("X1", "probability_1", colnames(pre_probability) )
colnames(pre_probability) <- gsub("X0", "probability_0", colnames(pre_probability) )
#dim(pre_probability) #[1] 1009 2
#colnames(data_test)
data_1 <- cbind(pre_probability, data_test[,c(ncol(data_test),ncol(data_test)-1, ncol(data_test)-2)])
head(data_1)
## probability_0 probability_1 Freq Virus Drug
## 1 1.000000e+00 1.791792e-10 0 HIV-2 arbidol
## 2 5.594001e-04 9.994406e-01 0 hCoV-19-Alpha propagermanium
## 3 2.347075e-19 1.000000e+00 0 hCoV-19-Delta telaprevir
## 4 1.000000e+00 5.566425e-55 0 HPV-11 tenofovir-alafenamide
## 5 1.000000e+00 1.659949e-15 0 HBV cobicistat
## 6 2.146014e-22 1.000000e+00 0 HCV-7 chloroquine
#dim(data_1)
#table(cut(pre_probability$probability_1,breaks = c(-0.1,.5,1)))
#head(pre_probability)
#write.table(pre_probability,paste0(dir_path,Sys.Date(),"-ecfp4_auc_nb_original.txt"),row.names = F,sep = "\t")
#write.csv(pre_probability,paste0(dir_path,Sys.Date(),"-ecfp4_auc_nb_original_0_52.csv"),row.names = F)
pred_auc <- prediction(data_1$probability_1,data_1$Freq)
auc_roc <- performance(pred_auc,"auc")@y.values[[1]]
auc_roc #[1] 0.7839947
## [1] 0.7839947
write.csv(data_1,paste0(dir_path,Sys.Date(),"-nb_upsample_probability.csv"),row.names = F)
#######################plot
#View(data_1)
head(data_1)
## probability_0 probability_1 Freq Virus Drug
## 1 1.000000e+00 1.791792e-10 0 HIV-2 arbidol
## 2 5.594001e-04 9.994406e-01 0 hCoV-19-Alpha propagermanium
## 3 2.347075e-19 1.000000e+00 0 hCoV-19-Delta telaprevir
## 4 1.000000e+00 5.566425e-55 0 HPV-11 tenofovir-alafenamide
## 5 1.000000e+00 1.659949e-15 0 HBV cobicistat
## 6 2.146014e-22 1.000000e+00 0 HCV-7 chloroquine
library(ggplot2)
data_1$Freq <- factor(data_1$Freq)
ggplot(data_1) +
stat_density(aes(x = probability_1, color = Freq),size = 1,
alpha=0.5, bw = 0.01, geom="line",position="identity") +
scale_y_discrete(expand = c(0.01, 0)) +
scale_x_continuous(expand = c(0, 0),limits = c(-0.05, 1.05),
breaks = seq(0,1,.2)) +
labs(colour="Antiviral drug",x = "Probability", y = "Density",title = "NB")+
geom_vline(aes(xintercept=0.5),
color="gray80", linetype="dashed", size=0.5)+
theme(panel.spacing = unit(0.1, "cm"),
legend.position= "top",
legend.key = element_rect(colour = NA, fill = NA),
legend.text=element_text(size=14),
legend.title = element_text(size=14),
axis.ticks = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.line = element_line(colour = "black",
size = 0.5, linetype = "solid"),
axis.text =element_text(face="plain", color="black", family = "sans",
size=14,angle = 0),
panel.background = element_rect(fill = "white",
colour = "white",
size = 0.5, linetype = "solid"),
panel.grid.major = element_line(size = 1, linetype = 'dashed',
colour = "white"),
axis.title = element_text(color="black", size=14, face="plain",family="sans")) +
scale_color_manual(name = "Antiviral drug",values=c("#3b58a7","#90278e"))
##########################output
ggsave(filename = paste0(Sys.Date(),"-probability_plot_SVM.tif"), plot = last_plot(),
device = "tiff", path = dir_path,
scale = 1, width = 16, height = 12, units = "cm",
dpi = 300, limitsize = TRUE, compression = "lzw")
#View(data_1)
data_11 <- data_1[data_1$Freq == 1, ]
data_22 <- data_1[data_1$Freq != 1, ]
#dim(data_1) #[1] 115 5
table(cut(data_11$probability_1, breaks = c(0,0.5,1), include.lowest = T))
##
## [0,0.5] (0.5,1]
## 14 101
table(cut(data_22$probability_1, breaks = c(0,0.5,1), include.lowest = T))
##
## [0,0.5] (0.5,1]
## 530 364