发帖
雷达卡
最近在做的R语言与数据分析课程论文，在此给出代码，使用数据集可在UCI机器学习数据库中查找（红葡萄酒），代码上传后会乱序，具体可见附件中的记事本文件。
##############################读取数据，缺失值处理##############################
library(caret)
dat0=read.csv('C:/Users/DELL/Desktop/redwines.csv',header = T)
summary(dat0)
nrow(dat0)
rw=na.omit(dat0)
nrow(rw)
# 统计字符变量占比
table(rw$quality)
prop.table(table(rw$quality))
#####################计算红葡萄酒各项指标的基本统计量##########################
outline=matrix(0,ncol=ncol(rw),nrow=13)
for(i in 1:ncol(rw))
{
x<- rw[,i]#将第i个指标的数据赋值给x
n<- length(x) #N为样本的个数
m<- mean(x) #Mean为样本的均值
v<- var(x) #Var为样本的方差
s<- sd(x) #std_dev为样本标准差
me<- median(x) #Median为样本中位数
cv<- 100*s/m #CV为样本的变异系数
css<- sum((x-m)^2) #CSS为样本校正平方和
uss<- sum(x^2) #USS为样本未校正平方和
R<- max(x)-min(x) #R为样本极差
R1<- quantile(x,3/4)-quantile(x,1/4) #R1为样本半极差
sm<- s/sqrt(n) #std_mean为样本的标准误
g1<- n/((n-1)*(n-2))*sum((x-m)^3)/s^3 #Skewness为样本峰度系数
g2<- ((n*(n+1))/((n-1)*(n-2)*(n-3))*sum((x-m)^4)/s^4- (3*(n-1)^2)/((n-2)*(n-3))) #Kurtosis为样本偏度系数
#将各项数值依次赋值给outline矩阵元素
outline[1,i]<-n
outline[2,i]<-m
outline[3,i]<-v
outline[4,i]<-s
outline[5,i]<-me
outline[6,i]<-cv
outline[7,i]<-css
outline[8,i]<-uss
outline[9,i]<-R
outline[10,i]<-R1
outline[11,i]<-sm
outline[12,i]<-g1
outline[13,i]<-g2
}
write.csv(outline,file="D:/outline1.csv") #将结果输出
########################探索红酒质量是否服从正态分布######################
#绘制直方图
hist(rw$quality,xlab='红葡萄酒质量',ylab='频数',main='红葡萄酒质量直方图与折线图')
x=c(3.25,3.75,4.75,5.75,6.75,7.75)
y=c(10,53,681,638,199,18)
#绘制红葡萄酒质量分布折线图
lines(x,y,type = "l")
#绘制正态QQ图
qqnorm(y)
qqline(y)
#正态性W检验
shapiro.test(y)
#################################数据分割#################################
#设置随机种子
set.seed(1234)
#将数据集的80%划分为训练集，20%划分为测试集
#createDataPartition会自动从y的各个level随机取出等比例的数据组成训练集
trainIndex=createDataPartition(rw$quality,p=.8, list=FALSE,times=1)
rw_train=rw[trainIndex,] #训练集
rw_test=rw[-trainIndex,] #测试集
#设定分类标签
rw_train_labels <-rw_train[, 12]
rw_test_labels <-rw_test[, 12]
##################################KNN分类###############################
library(class)
knn_pred <- knn(train = rw_train, test = rw_test,
cl = rw_train_labels, k=1)
#评估模型性能
library(gmodels)
#使用混淆矩阵查看模型分类结果
CrossTable(x = rw_test_labels, y = knn_pred,
prop.chisq = F, prop.c = F,prop.r = T,
dnn = c('ACTUAL RW','PREDICTED RW'))
# 将预测目标转换为数值型，才能计算auc
test_numeric <- as.numeric(rw_test_labels)
knn_pred_numeric <- as.numeric(knn_pred)
# 调用pROC包的roc函数
library(pROC)
knn_auc <- roc(test_numeric,knn_pred_numeric)
# 输出auc曲线面积
print(knn_auc)
# 画auc图
plot(knn_auc,ylim=c(0,1),print.thres=TRUE,main=paste('AUC',round(knn_auc$auc[[1]],2)))
abline(h=1,col="blue",lwd=2)
abline(h=0,col="red",lwd=2)
#提升模型质量，使用Z-score得分进行标准化
preProcValues=preProcess(rw_train, method = "scale")
rw_train_z=predict(preProcValues, rw_train)
rw_test_z=predict(preProcValues, rw_test)
#使用knn算法对标准化数据进行分类
knn_pred1 <- knn(train = rw_train_z, test = rw_test_z,
cl = rw_train_labels, k=1)
CrossTable(x = rw_test_labels, y = knn_pred1,
prop.chisq=F, prop.c = F,prop.r = T,
dnn = c('ACTUAL RW','PREDICTED RW'))
#绘制ROC曲线
knn_pred1_numeric <- as.numeric(knn_pred1)
knn_auc1 <- roc(test_numeric,knn_pred1_numeric)
print(knn_auc1)
plot(knn_auc1,ylim=c(0,1),print.thres=TRUE,main=paste('AUC',round(knn_auc1$auc[[1]],2)))
abline(h=1,col="blue",lwd=2)
abline(h=0,col="red",lwd=2)
#尝试不同的k值
#k=2
knn_pred2 <- knn(train = rw_train_z, test = rw_test_z, cl = rw_train_labels, k=2)
CrossTable(x = rw_test_labels, y = knn_pred2, prop.chisq=F, prop.c = F)
#k=3
knn_pred3 <- knn(train = rw_train_z, test = rw_test_z, cl = rw_train_labels, k=3)
CrossTable(x = rw_test_labels, y = knn_pred3, prop.chisq=F, prop.c = F)
#k=5
knn_pred4 <- knn(train = rw_train_z, test = rw_test_z, cl = rw_train_labels, k=5)
CrossTable(x = rw_test_labels, y = knn_pred4, prop.chisq=F, prop.c = F)
#k=7
knn_pred5 <- knn(train = rw_train_z, test = rw_test_z, cl = rw_train_labels, k=7)
CrossTable(x = rw_test_labels, y = knn_pred5, prop.chisq=F, prop.c = F)
#k=9
knn_pred6 <- knn(train = rw_train_z, test = rw_test_z, cl = rw_train_labels, k=9)
CrossTable(x = rw_test_labels, y = knn_pred6, prop.chisq=F, prop.c = F)
#################################朴素贝叶斯算法###############################
#在使用接下来的模型时，需要将结果变量转化为因子型。重新对数据进行整理
x=rw_train[,-12]
y=as.factor(rw_train$quality)
dat=data.frame(x,y)
x1=rw_test[,-12]
y1=as.factor(rw_test$quality)
dat1=data.frame(x1,y1)
#使用朴素贝叶斯算法进行分类
library(e1071)
Bay<- naiveBayes(x,y)
Bay_pred<-predict(Bay,x1)
CrossTable(rw_test_labels,Bay_pred, prop.chisq=F, prop.c = F)
#绘制ROC曲线
test_numeric <- as.numeric(rw_test_labels)
Bay_pred_numeric <- as.numeric(Bay_pred)
Bay_auc <- roc(test_numeric,Bay_pred_numeric)
print(Bay_auc)
plot(Bay_auc,ylim=c(0,1),print.thres=TRUE,main=paste('AUC',round(Bay_auc$auc[[1]],2)))
abline(h=1,col="blue",lwd=2)
abline(h=0,col="red",lwd=2)
#提升模型性能，通过修改laplace的数值,修改数值结果变动不大，此处只显示修改一次
Bay1<- naiveBayes(x,y,laplace = 1)
Bay1_pred<-predict(Bay1,x1)
Bay1_agree <- Bay1_pred == rw_test$quality
table(Bay1_agree)
prop.table(table(Bay1_agree))
######################################决策树##################################
set.seed(1234)
#C50决策树算法
library(C50)
C5<- C5.0(y ~ ., data = dat)
C5_pred<-predict(C5,dat1)
summary(C5)
#回归树
library(rpart)
rpart<-rpart(y~.,data=dat)
rpart
#可视化决策树
library(rpart.plot)
rpart.plot(rpart,digits=3)
rpart.plot(rpart,digits=4,fallen.leaves=TRUE,type=3,extra=101)
#C50决策树分类混淆矩阵
CrossTable(rw_test_labels,C5_pred, prop.chisq=F, prop.c = F)
#绘制ROC曲线图
test_numeric <- as.numeric(rw_test_labels)
C5_pred_numeric <- as.numeric(C5_pred)
C5_auc <- roc(test_numeric,C5_pred_numeric)
print(C5_auc)
plot(C5_auc,ylim=c(0,1),print.thres=TRUE,main=paste('AUC',round(C5_auc$auc[[1]],2)))
abline(h=1,col="blue",lwd=2)
abline(h=0,col="red",lwd=2)
#加入自适应增强算法
C5_boost100<-C5.0(y~.,data=dat,trials=100)
summary(C5_boost100)
C5_boost100_pred<-predict(C5_boost100,dat1)
#输出混淆矩阵
CrossTable(rw_test_labels,C5_boost100_pred, prop.chisq=F, prop.c = F)
#绘制ROC曲线图
test_numeric <- as.numeric(rw_test_labels)
C5_boost100_pred_numeric <- as.numeric(C5_boost100_pred)
C5_boost100_auc <- roc(test_numeric,C5_boost100_pred_numeric)
print(C5_boost100_auc)
plot(C5_boost100_auc,ylim=c(0,1),print.thres=TRUE,main=paste('AUC',round(C5_boost100_auc$auc[[1]],2)))
abline(h=1,col="blue",lwd=2)
abline(h=0,col="red",lwd=2)
################################支持向量机#####################################
#使用支持向量机进行分类
library(kernlab)
svm_van <- ksvm(y ~ ., data = dat,kernel = "vanilladot")
#查看模型基础信息
svm_van
#利用模型进行预测
svm_pred <- predict(svm_van, dat1)
head(svm_pred)
#将测试集中的预测数据与真实数据进行比较
table(svm_pred, rw_test$quality)
agree <- svm_pred == rw_test$quality
table(agree)
prop.table(table(agree))
#提高模型性能
#径向基“高斯”核
svm_rbf <- ksvm(y ~ ., data = dat, kernel = "rbfdot")
svm_rbf
svm_rbf_pred <- predict(svm_rbf, dat1)
svm_rbf_agree <- svm_rbf_pred== rw_test$quality
table(svm_rbf_agree)
prop.table(table(svm_rbf_agree))
#多项式内核
svm_poly<- ksvm(y ~ ., data = dat, kernel = "polydot")
svm_poly
svm_poly_pred <- predict(svm_poly, dat1)
svm_poly_agree <- svm_poly_pred == rw_test$quality
table(svm_poly_agree)
prop.table(table(svm_poly_agree))
#双曲正切核
svm_tan<- ksvm(y ~ ., data = dat, kernel = "tanhdot")
svm_tan
svm_tan_pred <- predict(svm_tan, dat1)
svm_tan_agree <- svm_tan_pred == rw_test$quality
table(svm_tan_agree)
prop.table(table(svm_tan_agree))
#拉普拉斯核
svm_lap<- ksvm(y ~ ., data = dat, kernel = "laplacedot")
svm_lap
svm_lap_pred <- predict(svm_lap, dat1)
svm_lap_agree <- svm_lap_pred == rw_test$quality
table(svm_lap_agree)
prop.table(table(svm_lap_agree))
#贝塞尔内核
svm_bes<- ksvm(y ~ ., data = dat, kernel = "besseldot")
svm_bes
svm_bes_pred <- predict(svm_bes, dat1)
svm_bes_agree <- svm_bes_pred == rw_test$quality
table(svm_bes_agree)
prop.table(table(svm_bes_agree))
#ANOVA RBF内核
svm_ano<- ksvm(y ~ ., data = dat, kernel = "anovadot")
svm_ano
svm_ano_pred <- predict(svm_ano, dat1)
svm_ano_agree <- svm_ano_pred == rw_test$quality
table(svm_ano_agree)
prop.table(table(svm_ano_agree))
#样条内核
svm_spl<- ksvm(y ~ ., data = dat, kernel = "splinedot")
svm_spl
svm_spl_pred <- predict(svm_spl, dat1)
svm_spl_agree <- svm_spl_pred == rw_test$quality
table(svm_spl_agree)
prop.table(table(svm_spl_agree))
#使用高斯核函数进行预测
CrossTable(x = rw_test_labels, y = svm_rbf_pred,
prop.chisq = F, prop.c = F,prop.r = T,
dnn = c('ACTUAL RW','PREDICTED RW'))
#绘制ROC曲线图
test_numeric <- as.numeric(rw_test_labels)
svm_rbf_pred_numeric <- as.numeric(svm_rbf_pred)
svm_rbf_auc <- roc(test_numeric,svm_rbf_pred_numeric)
print(svm_rbf_auc)
plot(svm_rbf_auc,ylim=c(0,1),print.thres=TRUE,main=paste('AUC',round(svm_rbf_auc$auc[[1]],2)))
abline(h=1,col="blue",lwd=2)
abline(h=0,col="red",lwd=2)
####################################随机森林###################################
library(randomForest)
set.seed(1234)
rf<- randomForest(y ~.,data=dat,ntree=100,
mtry=sqrt(ncol(dat)), maxnotes=10)
rf
# type可以是"response","prob","vote",分别表示输出预测向量是预测类别、预测概率或投票矩阵
rf_pred <- predict(rf, dat1, type='response')
rf_agree <- rf_pred == rw_test$quality
table(rf_agree)
prop.table(table(rf_agree))
CrossTable(x = rw_test_labels, y = rf_pred,
prop.chisq = F, prop.c = F,prop.r = T,
dnn = c('ACTUAL RW','PREDICTED RW'))
#绘制ROC曲线图
test_numeric <- as.numeric(rw_test_labels)
rf_pred_numeric <- as.numeric(rf_pred)
rf_auc <- roc(test_numeric,rf_pred_numeric)
print(rf_auc)
plot(rf_auc,ylim=c(0,1),print.thres=TRUE,main=paste('AUC',round(rf_auc$auc[[1]],2)))
abline(h=1,col="blue",lwd=2)
abline(h=0,col="red",lwd=2)
京公网安备 11010802022788号
