#############R分类
###什么是分类问题
###创建误差|分类混淆矩阵
cp <- read.csv("college-perf.csv")
str(cp)
#因子重新排序
dput(levels(cp$Perf))
#常规因子排序
cp$Perf <- factor(cp$Perf,levels=c("Low", "Medium","High"))
str(cp)
#顺序因子
cp$Pred <- ordered(cp$Pred,levels=c("Low", "Medium","High"))
cp$Perf <- ordered(cp$Perf,levels=c("Low", "Medium","High"))
str(cp)
#reorder()排序
# 使用函数reorder()，该函数有三个参数：因子、排序依据的数据和汇总数据的函数
# 复制一份数据
iss <- InsectSprays
iss$spray
iss$spray <- reorder(iss$spray,iss$count,FUN = mean)
iss$spray
#混淆矩阵
tab <- table(cp[,c('Perf','Pred')])
tab
#百分比
prop.table(tab)
round(prop.table(tab,1)*100,1)
#将向量转为数据框
a <- paste(round(prop.table(tab,1)*100,1),'%',sep='')
paste0(round(prop.table(tab,1)*100,1),'%')
as.data.frame(matrix(a,ncol=3,dimnames=list(Perf=c("Low", "Medium","High"),Pred=c("Low", "Medium","High"))))
#混淆矩阵可视化
barplot(tab,legend=TRUE)
mosaicplot(tab,main='prediction performance')

#在不同类别上比较模型的表现
summary(tab)

library(caret)
confusionMatrix(tab)
?confusionMatrix


###创建ROC图
library(ROCR)
dat <- read.csv("roc-example-1.csv")
head(dat)
#创建预测对象
pred <- prediction(dat$prob,dat$class)
#创建预测表现对象
?performance
perf <- performance(pred,'tpr','fpr')
#画图
old.par <- par()
par(mfrow=c(1,2))
plot(perf)
lines(par()$usr[1:2],par()$usr[3:4])
#利用基准值
perf2 <- performance(pred,'acc','cutoff')
str(perf2)
plot(perf2)
par(old.par)

#查看AUC值
performance(pred,'auc')

#找出各个真阳性比率的截断点
str(perf)
prob.cuts <- data.frame([email protected][[1]],[email protected][[1]],[email protected][[1]])
head(prob.cuts)
tail(prob.cuts)
#选择真阳性比率对应的截断值


#使用任意的分类标签
dat <- read.csv("roc-example-2.csv")
head(dat)
pred <- prediction(dat$prob,dat$class,label.ordering = c('non-buyer','buyer'))#通过传递分类标签分量明确指出标签
perf <- performance(pred,'tpr','fpr')
plot(perf)
lines(par()$usr[1:2],par()$usr[3:4])


###构建、绘制、和评估——分类树
library(rpart)
library(rpart.plot)
library(caret)
bn <- read.csv("banknote-authentication.csv")
head(bn)
str(bn)

#创建数据分块
set.seed(1000)
table(bn$class)
train.idx <- createDataPartition(bn$class,p=0.7,list=FALSE)

#创建树
?rpart
mod <- rpart(class~.,data=bn[train.idx,],method='class',control=rpart.control(minsplit=20,cp=0.01))

#查看文本输出
mod

#创建树图
prp(mod,type=4,extra=2,digits=3)
prp(mod,type=2,extra=104,nn=TRUE,fallen.leaves = TRUE,faclen=4,varlen=8,shadow.col = 'gray')
#faclen分类名
#varlen变量名

plot(mod,compress=TRUE,margin=0.2)
text(mod,cex=0.6)

#给树剪枝,防止过度拟合
str(mod)
mod$cptable #xerror is not greater than  minimum xerror + xstd
mod.pruned <- prune(mod,cp=mod$cptable[6,'CP'])#截断复杂性水平
0.1276596+0.01687712

#查看剪枝后的树
prp(mod.pruned,type=2,extra=104,nn=TRUE,fallen.leaves=TRUE,faclen=4,varlen=8,shadow.col='gray')

#用剪枝后的模型对验证模块做预测
pred.pruned <- predict(mod.pruned,bn[-train.idx,],type='class')

#创建误差/分类-混淆矩阵:可视化
table(bn[-train.idx,]$class,pred.pruned,dnn=c('actual','predicted'))

#计算原始概率
pred.pruned <- predict(mod.pruned,bn[-train.idx,],type='prob')
#创建ROC图
library(ROCR)
pred <- prediction(pred.pruned[,2],bn[-train.idx,]$class)
perf <- performance(pred,'tpr','fpr')
plot(perf)
abline(a=0,b=1)
lines(par()$usr[1:2],par()$usr[3:4])

#补充知识：https://blog.csdn.net/yujunbeta/article/details/12648717
# 1、CP参数定义
# 
# cp: complexity parameter复杂性参数，用来修剪树的.
# 
# 当决策树复杂度超过一定程度后，随着复杂度的提高，测试集的分类精确度反而会降低。因此，建立的决策树不宜太复杂，
# 需进行剪枝。该剪枝算法依赖于复杂性参数cp,cp随树复杂度的增加而减小，当增加一个节点引起的分类精确度变化量小于树
# 复杂度变化的cp倍时，则须剪去该节点。故建立一棵既能精确分类，又不过度适合的决策树的关键是求解一个合适的cp值


# 2、CP选择问题。
# 
# 建立树模型要权衡两方面问题，一个是要拟合得使分组后的变异较小，另一个是要防止过度拟合，而使模型的误差过大，
# 前者的参数是CP，后者的参数是Xerror。所以要在Xerror最小的情况下，也使CP尽量小。

###条件推断决策树party::ctree



###随机森林:randomForest 
library(randomForest)
library(caret)
bn <- read.csv("banknote-authentication.csv")
head(bn)
##将响应变量转为因子
bn$class <- factor(bn$class)
str(bn)
?randomForest

mod <- randomForest(x=bn[,1:4],y=bn[,5],ntree=500)
mod
pred <-predict(mod,bn)
table(bn$class,pred,dnn=c('acutal','predicted'))

##参数调优
#方法之一交叉检验
#主要参数：
#决策树棵数ntree
#确定节点分隔标准时要考虑的变量个数mtry
#创建组合列表：expand.grid()，返回值组合所有水平的数据框，参数可以是向量、因子或保存其的列表
grid <- expand.grid(ntree=c(10,100,200),mtry=c(3,4))
#评估模型在各种组合参数下的性能，重复执行3次
library(cvTools)
library(foreach)
set.seed(719)
K = 10
R = 3
cv <- cvFolds(NROW(bn),K=K,R=R)
grid <- expand.grid(ntree=c(10,100,1000),mtry=c(3,4))

result <- foreach(g=1:NROW(grid),.combine=rbind) %do% {
  foreach(r=1:R,.combine = rbind) %do% {
    foreach(k=1:K,.combine=rbind) %do%{
      val.idx <- cv$subsets[which(cv$which==k),r]
      train <-bn[-val.idx,]
      test <- bn[val.idx,]
      #模型训练
      m <- randomForest(class~.,data=train,ntree=grid[g,'ntree'],mtry=grid[g,'mtry'])
      #预测
      predicted <- predict(m,newdata=test)
      #性能评估
      precision <- sum(predicted==test$class)/NROW(predicted)
      return(data.frame(g=g,precision=precision))
    }
  }
}
result
#查看每个g值绑定的平均性能
library(plyr)
ddply(result,.(g),summarize,mean_precision=mean(precision))
#查找高性能的g与之对应的参数组合
grid[2,]

m<-randomForest(x=bn[,1:4],y=bn[,5],ntree=100,mtry=3,keep.forest=TRUE)
str(m)

bn <- read.csv("banknote-authentication.csv")
bn$class <- factor(bn$class)
set.seed(1000)
sub.idx <- createDataPartition(bn$class, p=0.7, list=FALSE)
mod <- randomForest(x = bn[sub.idx,1:4], y=bn[sub.idx,5],ntree=500, keep.forest=TRUE)
pred <- predict(mod, bn[-sub.idx,])
table(bn[-sub.idx,"class"], pred, dnn = c("Actual", "Predicted"))
#计算原始概率，创建ROC图
probs <- predict(mod,bn[-sub.idx,],type='prob')
#创建ROC曲线图
pred <- prediction(probs[,2],bn[-sub.idx,'class'])
perf <- performance(pred,'fpr','tpr')
plot(perf)
#指定分类截断值



##评估变量的重要性
#正则化随机森林 RRF包
m <- randomForest(class~.,data=bn,importance=TRUE)
importance(m)#准确度与基尼不纯度
varImpPlot(m)

install.packages('RRF')
library(RRF)
??RRF

###用支持向量机分类
###支持向量机分类
library(e1071)
library(caret)
bn <- read.csv("banknote-authentication.csv")
str(bn)
#将输入变量转为因子
bn$class <- factor(bn$class)
#为数据分块
set.seed(1000)
t.idx <- createDataPartition(bn$class,p=0.7,list=FALSE)
#建立模型
?svm
mod <- svm(class~.,data=bn[t.idx,])
mod
#建立误差/分类——混淆矩阵验证模型训练集上的表现
table(fitted(mod),bn[t.idx,'class'],dnn=c('pred','acu'))
#在验证模型上考察模型的表现
table(predict(mod,bn[-t.idx,]),bn[-t.idx,'class'],dnn=c('pre','acu'))
#在训练分上绘制模型的图形#选择变量
plot(mod,data=bn[t.idx,],skew ~ variance)
#在验证分块上绘制图形
plot(mod,data=bn[-t.idx,],skew ~ variance)

#控制变量缩放:所有变量（x,y)
scale()
#决定SVM模型的类型:type设置

#分类别赋予权重class.weights=c('0'=0.3,'1'=0.7)
table(bn$class)
mod <- svm(class~.,data=bn[t.idx,],class.weights=c('0'=0.3,'1'=0.7))

##参数调优
#方式一：交叉验证
#方式二：使用包提供函数
?tune()
tune.svm(Species~., data = iris, gamma = 2^(-1:1), cost = 2^(2:4))
svm(x,y,scale=TRUE,type=NULL,kernel='radial',gamma=1,cost=1)



###用朴素贝叶斯分类
library(e1071)
library(caret)
ep <- read.csv("electronics-purchase.csv")
head(ep)
str(ep)
#要求所有变量为分类变量
#为数据分块
set.seed(1000)
t.idx <- createDataPartition(ep$Purchase,p=0.67,list=FALSE)
#建立模型
?naiveBayes
mod <- naiveBayes(Purchase~.,data=ep[t.idx,])
mod
#在验证模型上考察模型的表现
table(predict(mod,ep[-t.idx,]),ep[-t.idx,'Purchase'],dnn=c('pre','acu')


###用k近邻分类
library(class)
library(caret)
#要求自变量为数值类型，应变量为分类变量
vac <- read.csv("vacation-trip-classification.csv")
head(vac)
str(vac)
#将预测变量标准化
vac$Income <- scale(vac$Income)
vac$Family_size <- scale(vac$Family_size)
#将数据分块
set.seed(1000)
train.idx <- createDataPartition(vac$Result, p = 0.5, list = FALSE)
train <- vac[train.idx, ]
temp <- vac[-train.idx, ]
val.idx <- createDataPartition(temp$Result, p = 0.5, list = FALSE)
val <- temp[val.idx, ]
test <- temp[-val.idx, ]
#k=1时，生成验证集的结果
predicted<-knn(train[,1:2],val[,1:2],train[,3],k=1)
#对k=1生成误差矩阵
table(val$Result,predicted,dnn=c('acu','pre'))
#重复多种K
knn.automate <- function(train.var,val.var,train.fac,val.trg,start_k,end_k) {
  for(k in start_k:end_k) {
    pred <- knn(train.var,val.var,train.fac,k)
    tab <- table(val.trg,pred,dnn=c('acu','pred'))
    cat('error matrix for k = ',k,'\n')
    cat("===============\n")
    print(tab)
    cat('------------\n\n\n')
  }
}

knn.automate(train[,1:2],val[,1:2],train[,3],val[,3],1,7)

#最佳K
#最佳k对预测分块中样本预测
pred.test <-  knn(train[,1:2],test[,1:2],train[,3],k=4)
table(test$Result,pred.test,dnn=c('acu','pred'))

#用knn计算原始概率而非预测值
#有助类失衡数据，计算原始概率选择一个用于分类的阶段概率
knn(train[,1:2],val[,1:2],train[,3],k=4,prob=TRUE)

###用神经网络分类
library(nnet)
library(caret)
bn <- read.csv("banknote-authentication.csv")
#要求所有自变量为数值类型，应变量为0-1型
#将输出变量转为因子
bn$class <- factor(bn$class)
#为数据分块
tra.idx <- createDataPartition(bn$class,p=0.7,list=FALSE)
#建立模型
mod <- nnet(class~.,data=bn[tra.idx,],size=3,maxit=10000,decay=0.001,rang=0.05)
#size:初略估计隐匿层单元数量设置为输入层单元数量和输出层单元数量的均值，更大的值会得到更好的结果，但代价是计算时间变得更长
#maxit:算法收敛而做的最大迭代次数
#decay:权重衰减
#rang:指定初始随机权重范围[-rang,rang],如果输入值很大，选择能令rang*(max|variable|)接近于1的权重范围。
###更好控制nnet的其他选项
#na.action=na.omit:来排除包含缺失值的样本
#skip=TRUE:添加从输入节点到输出节点的直接连接
#linout=FALSE:默认值为FALSE,输出层不使用线性函数而使用sigmoid函数
#entropy=FALSE:默认使用SSE.预测对象的分类水平为2时，使用信息熵推断参数，分类数大于3时，使用SSE推断参数，并使用softmax.
#softmax=FALSE:神经网络输出为概率形式，应使用softmax函数就爱那个最小值变为0，最大值变为1.
#MaxNWts=1000:默认加权值的最大个数为1000。 

#用模型来预测验证分块
pred <- predict(mod,newdata=bn[-tra.idx,],type='raw') #class表示预测值为分类，raw表示预测值为原始概率

##直接指定X与Y
class.ind(bn$class)
mod <- nnet(bn[tra.idx,1:4],class.ind(bn[tra.idx,'class']),size=3,entropy=TRUE)
pred <- predict(mod,newdata=bn[-tra.idx,],type='class')
table(bn[-tra.idx,'class'],pred)


###用逻辑回归分类
library(caret)
bh <- read.csv("boston-housing-logistic.csv")
#要求所有的自变量都是数值型的,因变量为二分类变量
#将输出变量转化为因子:glm函数会默认把第一个因子水平看作失败,将其他水平视为成功
bh$CLASS <- factor(bh$CLASS,levels=c('0','1'))
str(bh)
#为数据分块
set.seed(1000)
tra.idx <- createDataPartition(bh$CLASS,p=0.7,list=FALSE)
#建立回归模型
logit <- glm(CLASS~.,data=bh[tra.idx,],family='binomial')
#查看模型
summary(logit)
#deciance residuals为对数成败比率而非概率值的变异分布
#计算验证分块中样本的成功概率并保存变量
bh[-tra.idx,'prob_succ'] <- predict(logit,bh[-tra.idx,],type='response') 
#用0.5的阶段概率将样本分类
bh[-tra.idx,'prob_50'] <- ifelse(bh[-tra.idx,'prob_succ']<0.5,0,1)
#建立分类混淆矩阵
table(bh[-tra.idx,'prob_50'],bh[-tra.idx,'CLASS'],dnn=c('pre',"acu"))

Recipe: Classifying using Linear Discriminant Function Analysis
#--------------------------------------------------------------
##用线性判别函数分类
library(MASS)
library(caret)

bn <- read.csv("banknote-authentication.csv")

bn$class <- factor(bn$class)

set.seed(1000)
t.idx <- createDataPartition(bn$class, p = 0.7, list=FALSE)

ldamod <- lda(bn[t.idx, 1:4], bn[t.idx, 5])

bn[t.idx,"Pred"] <- predict(ldamod, bn[t.idx, 1:4])$class

table(bn[t.idx, "class"], bn[t.idx, "Pred"], dnn = c("Actual", "Predicted"))

bn[-t.idx,"Pred"] <- predict(ldamod, bn[-t.idx, 1:4])$class 

table(bn[-t.idx, "class"], bn[-t.idx, "Pred"], dnn = c("Actual", "Predicted"))

ldamod <- lda(class ~ ., data = bn[t.idx,])



#Recipe: Using AdaBoost to combine classification tree models
#---------------------------
##用AdaBoost来整合分类树模型
library(caret)
library(ada)

bh <- read.csv("banknote-authentication.csv")

bn$class <- factor(bn$class)
set.seed(1000)
t.idx <- createDataPartition(bn$class, p=0.7, list=FALSE)

cont <- rpart.control()
#创建默认对象

mod <- ada(class ~ ., data = bn[t.idx,], iter=50, loss="e", type="discrete", control = cont)
#ada会使用rpart函数生成多棵分类数，需要提供rpart.control对象。
#iter=50;设置增强迭代的上限
#loss='e',设置损失函数为指数函数
#type=discrete :指定模型类型为分类非回归

mod

pred <- predict(mod, newdata = bn[-t.idx,], type = "vector")

table(bn[-t.idx, "class"], pred, dnn = c("Actual", "Predicted"))