rm(list = ls()); ls()
## character(0)
library(nnet)
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
library(mltools)
###############################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
###############################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
#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)
## 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
#str(train_data$num)
##model
#####################nnet_parameter
#newdata_tune_nnet <- tune.nnet(num ~ .,data = train_data, size = 1:3,
# decay = c(0.0001,0.001,0.01,0.1),maxit = c(10,100,200,500,800,1000),
# tunecontrol = tune.control(nrepeat = 5))
#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 = "over")$data
fit_nnet_ori <- nnet(num ~ .,data = traindata_ori,size=10,MaxNWts = 100000)
## # weights: 1491
## initial value 1982.283082
## iter 10 value 1403.617627
## iter 20 value 618.384443
## iter 30 value 223.492690
## iter 40 value 123.305847
## iter 50 value 86.610442
## iter 60 value 63.088237
## iter 70 value 51.682517
## iter 80 value 42.694403
## iter 90 value 22.904809
## iter 100 value 16.942807
## final value 16.942807
## stopped after 100 iterations
#########################################################
#### auc_rf
n_col_newdata <- ncol(data_train)
pre_auc_nnet_ori <- predict(fit_nnet_ori, testdata_ori[,-n_col_newdata], type = 'raw')
pred_auc_nnet_ori <- prediction(pre_auc_nnet_ori[,1],testdata_ori[,n_col_newdata])
auc_nnet_ori <- performance(pred_auc_nnet_ori,"auc")@y.values[[1]]
auc_nnet_ori #[1] 0.9793877
## [1] 0.9178512
#head(pre_probability_ori)
modelroc_nnet_ori <- pROC::roc(testdata_ori[,n_col_newdata],as.numeric(pre_auc_nnet_ori[,1]))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
cutoff_nnet_ori <- coords(modelroc_nnet_ori, "best", ret = "threshold") #[1] 0.265
predict.results_nnet_ori <- ifelse(pre_auc_nnet_ori[,1] > as.numeric(cutoff_nnet_ori[1]),"1","0")
freq_default_nnet_ori <- table(predict.results_nnet_ori,testdata_ori[,n_col_newdata])
Sensitivity_nnet_ori <- freq_default_nnet_ori[1,1]/sum(freq_default_nnet_ori[,1])
Specificity_nnet_ori <- freq_default_nnet_ori[2,2]/sum(freq_default_nnet_ori[,2])
accuracy_nnet_ori <- mean(testdata_ori[,n_col_newdata] == predict.results_nnet_ori)
Balanced_accuracy_nnet_ori <- (Sensitivity_nnet_ori+Specificity_nnet_ori)/2
Balanced_accuracy_nnet_ori #[1] 0.954963
## [1] 0.8985149
library(caret)
## 载入需要的程辑包：ggplot2
#confusionMatrix(factor(predict.results_nnet_ori),factor(testdata_ori[,n_col_newdata]))
preds_nnet <- as.numeric(as.character(predict.results_nnet_ori))
actuals_nnet <- as.numeric(as.character(testdata_ori$num))
mcc_result_nnet <- mcc(actuals_nnet,preds_nnet)
mcc_result_nnet
## [1] 0.7454343
#####0 1
#0 1350 15
#1 116 13
#################################################################################################
##nb
data_train <- ovun.sample(num ~ ., data = data_train, "over", seed = 1)$data
dim(data_train) #[1] 4129 148
## [1] 4041 148
#################################################probability
#fit_probability <- svm(num ~ .,data = data_train, scale = F, probability=TRUE )
dim(data_train) #[1] 4981 116
## [1] 4041 148
table(data_train$num) #4515 466
##
## 0 1
## 2057 1984
fit_probability <- nnet(num ~ .,data = data_train,size=10,MaxNWts = 100000)
## # weights: 1491
## initial value 3003.200405
## iter 10 value 1850.316174
## iter 20 value 627.729778
## iter 30 value 287.064420
## iter 40 value 178.225116
## iter 50 value 114.585927
## iter 60 value 76.443542
## iter 70 value 60.491964
## iter 80 value 48.042932
## iter 90 value 39.243748
## iter 100 value 33.490860
## final value 33.490860
## stopped after 100 iterations
#### 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) <- c("probability_1")
#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)
#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_nnet_original.txt"),row.names = F,sep = "\t")
#write.csv(pre_probability,paste0(dir_path,Sys.Date(),"-ecfp4_auc_nnet_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.9207519
## [1] 0.9207519
write.csv(data_1,paste0(dir_path,Sys.Date(),"-nnet_upsample_probability.csv"),row.names = F)
#################################################plot
#View(data_1)
head(data_1)
## probability_1 Freq Virus Drug
## 1 0 0 HIV-2 arbidol
## 2 0 0 hCoV-19-Alpha propagermanium
## 3 0 0 hCoV-19-Delta telaprevir
## 4 0 0 HPV-11 tenofovir-alafenamide
## 5 0 0 HBV cobicistat
## 6 0 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 = "NNET")+
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]
## 20 95
table(cut(data_22$probability_1, breaks = c(0,0.5,1), include.lowest = T))
##
## [0,0.5] (0.5,1]
## 872 22