##数据预处理
df<-read.csv(file="旅游板块和航空板块数据(20210210).csv")         
dim(df)
head(df)
X<-df[,-1]
y<-c(rep(0,25),rep(1,29))
n<-length(y)
rownames(X)<-df[,1]
group1<-1:20
group2<-26:50
X.train1<-X[group1,]
X.train2<-X[group2,]
X.train<-as.matrix(rbind(X.train1,X.train2))
y.train<-y[c(group1,group2)]
n1<-length(y.train)
X.test<-as.matrix(X[-c(group1,group2),])
y.test<-y[-c(group1,group2)]

##逻辑回归
logist.fit<-glm(y.train~.,data=data.frame(X.train),family=binomial(link="logit"))
summary(logist.fit)
pre.nei<-predict(logist.fit,type = "response")
pre.wai<-predict(logist.fit,newdata=data.frame(X.test),type = "response")

#敏感性和特异性分析
#install.packages("pROC")
library(pROC) 
roc.fit <- roc(y.train,pre.nei) 
plot(roc.fit, print.auc=TRUE, auc.polygon=TRUE,legacy.axes=TRUE, grid=c(0.1, 0.2), 
     grid.col=c("green", "red"), max.auc.polygon=TRUE, auc.polygon.col="skyblue", print.thres=TRUE,xlab="特异性",ylab="敏感性") 

#阈值的确定
sum(diag(table(ifelse(pre.nei>0.611,1,0),y.train)))/length(y.train)
table(ifelse(pre.wai>0.611,1,0),y.test)/length(y.test)
ifelse(pre.wai>0.611,1,0)

##LASSO变量选择
#install.packages("glmnet")
library(glmnet)
cv.lasso<-cv.glmnet(X.train,y.train,alpha = 1,intercept = TRUE,standardize=FALSE,family=binomial(link="logit"))
plot(cv.lasso,xlab="log(λ)",ylab="广义线性模型偏差")
#cv.lasso$lambda[which.min(cv.lasso$cvm)]
lasso.fit<-glmnet(X.train,y.train,alpha = 1,lambda=cv.lasso$lambda.1se,
                  intercept = TRUE,standardize=FALSE,family=binomial(link="logit"))
lasso.beta<-coef(lasso.fit)
lasso.beta
lasso.nei<-1/(1+exp(-cbind(1,X.train)%*%lasso.beta))
lasso.wai<-1/(1+exp(-cbind(1,X.test)%*%lasso.beta))

sum(diag(table(ifelse(lasso.nei>0.611,1,0),y.train)))/length(y.train)
table(ifelse(lasso.wai>0.611,1,0),y.test)/length(y.test)
ifelse(lasso.wai>0.611,1,0)

##支持向量机
#install.packages("e1071")
library(e1071)
#线性
y.train1<-as.factor(ifelse(y.train==0,-1,1))
cv.svm<-tune.svm(X.train,y.train1, type= "C-classification",
         scale=FALSE, kernel="linear", cost=c(0.001,0.01,0.1,1,5,10,100,1000))
summary(cv.svm)
svm.fit<-cv.svm$best.model
svm.fit<-svm(y.train1~.,data=data.frame(X.train), type= "C-classification",
                 scale=FALSE, kernel="linear", cost=1, gamma=0.1)
summary(svm.fit)

svm.nei<-predict(svm.fit)
svm.wai<-predict(svm.fit,X.test)


sum(diag(table(svm.nei,y.train)))/length(y.train)
table(svm.wai,y.test)/length(y.test)
svm.wai

#非线性
ncv.svm<-tune.svm(X.train, y.train1, type= "C-classification",scale=FALSE, kernel="radial",
                  gamma=c(0.1,0.5,1,2,3,4),cost=c(0.001,0.01,0.1,1,5,10,100,1000))
summary(ncv.svm)
plot(ncv.svm,xlab="参数",ylab="成本",main="")
nsvm.fit<-ncv.svm$best.model
summary(nsvm.fit)
nsvm.nei<-predict(nsvm.fit)
nsvm.wai<-predict(nsvm.fit,X.test)

sum(diag(table(nsvm.nei,y.train)))/length(y.train)
table(nsvm.wai,y.test)/length(y.test)
nsvm.wai