该数据集最初来自国家糖尿病/消化/肾脏疾病研究所。数据集的目标是基于数据集中包含的某些诊断测量来诊断性的预测患者是否患有糖尿病。这里的所有患者都是Pima印第安至少21岁的女性。数据集由多个医学预测变量和一个目标变量组成Outcome。预测变量包括患者的怀孕次数、BMI、胰岛素水平、年龄等。
library(mlbench)
data(PimaIndiansDiabetes2)
data<-PimaIndiansDiabetes2
attach(data)
head(data,3)
summary(data)
library(mice)
md.pattern(data)
因为数据有缺失值，需要对数据先进行探索，先用caret包中的预处
理函数preProcess完成缺失值的插补后，再利用分类算法建立预测模型
，识别糖尿病患者。
preProcess函数可以对特征变量施行很多操作，包括中心化和标准化。
center 中心化，即减去自变量的平均值
scale 标准化，即除以自变量的标准差，变换后值域为[0,1]。如果新
样本值大于或小于训练集中值，则值将超出此范围。
library(caret)
#标准化
preprovalue<-preProcess(data[,-9],method=c("center","scale"))
fitted是拟合值，predict是预测值。模型是基于给定样本的值建立的，在这些给定样本上做预测就是拟合。在新样本上做预测就是预测。
举例：
fit<-lm(weight~height,data=women)
fitted(fit)
predict(fit,newdata=data.frame(height=90))
scaleddata<-predict(preprovalue,data[,-9])
caret包中的可视化是对lattice包作图的集成。主要函数是featurePlo
t:即featurePlot(x, y, plot = “strip”, labels = c("Feature",
""), ...)
x 为一个连续数值的矩阵或data frame;y 是一个表明类别的因子变量
plot 是作图的种类。对于分类有： box, strip,
density,pairs，ellipse，对于回归有pairs，scatter。默认情况下，
当y为因子，则为strip，否则为scatter
featurePlot(scaleddata,data[,9],plot = "box")
YeoJohnson转换与BoxCox变换相似，但是它的自变量可以是0或负数，
而BoxCox只能是正数
preprobox<-preProcess(scaleddata,method=c("YeoJohnson"))
boxdata<-predict(preprobox,scaleddata)
利用袋装算法进行缺失值插补(只能对数值型变量处理)
preprocimp<-preProcess(boxdata,method="bagImpute")
procdata<-predict(preprocimp,boxdata)
procdata$class<-data[,9]
决策树模型的基本思想，先从n个自变量中挑选1个，寻找最佳分割点，
将数据划分为两组，针对分组后数据将上述步骤重复下去，直到满足某
种条件。
在R中通常可以用rpart包来实现CART算法，其中重要的参数为cp(复杂参数，越小越复杂）,cp全称为complexity parameter，指某个点的复杂度，对每一步拆分,模型的拟合优度必须提高的程度，由control控制，下面将前面处理过的数据输入决策树模型，下面的cp参数设置为0，是为了让模型变得复杂，为方便后面演示剪枝处理
library(rpart)
rpartmodel<-rpart(class~.,data=procdata,control=rpart.control(cp=0))
plotcp(rpartmodel)
printcp(rpartmodel)
根据上面的输出自动求出最小的CP值，再用prune函数对树模型进行修剪
prune函数可以实现最小代价复杂度剪枝法，对于CART的结果，每个节点均输出一个对应的cp
prune函数通过设置cp参数来对决策树进行修剪,cp为复杂度系数
cptable<-as.data.frame(rpartmodel$cptable)
cptable$errsd<-cptable$xerror+cptable$xstd
#选择具有最小xerror的cp
cpvalue<-cptable[which.min(cptable$errsd),"cp"]
#prunemodel存储剪枝后的模型
prunemodel<-prune(rpartmodel,cpvalue)
#使用rpart.plot包来画出树结构图
library(rpart.plot)
rpart.plot(prunemodel)
还可以根据建立的树模型判断各解释变量的重要性，varImp函数可以根据树模型输出重要性度量，从下面数据可以看到，glucose即体内葡萄糖是最重要的变量，其他的变量如insulin胰岛素和年龄和体重与患病与否也有很重要 关系。
varImp(prunemodel)
除了解释之外，我们还可以用该模型来计算预测值，并和本身的真实值进行比对，即混淆矩阵
pre<-predict(prunemodel,procdata,type="class")
pretable<-table(pre,procdata$class)
accuracy<-sum(diag(pretable))/sum(pretable)
上表纵轴是预测值，横轴是真实值，落在对角线上的数字为预测正确的样本，由此可计算模型的准确率为84%。
进一步可计算灵敏度和特异度，灵敏度即真实为阴性条件下预测正确的比率，特异度为真实为阳性条件下预测正确的比率
pretable[1,1]/sum(pretable[,1])#灵敏度
pretable[2,2]/sum(pretable[,2])#特异度
观察到一共有120例是错误预测的样本，这两类错判的意义不一样，55例是本来有病但未发现，而65例是无病但误诊为患病，这两类失误发现的成本可能是不一样的。
如果同样一个模型，对它既进行训练又进行预测，显然不合适，这样往往会高估模型的准确性，导致过度拟合，一般多重交叉验证或者对模型的复杂度进行约束，即添加惩罚项
num<-sample(1:10,nrow(procdata),replace=T)
res<-array(0,dim=c(2,2,10))
n<-ncol(procdata)
for(i in 1:10){
train<-procdata[num!=i,]
test<-procdata[num==i,]
model<-rpart(class~.,data = train,control=rpart.control(cp=0.1))
pre<-predict(model,test[,-n],type="class")
res[,,i]<-as.matrix(table(pre,test[ ,n]))
}
table<-apply(res,MARGIN = c(1,2),sum)
rate<-sum(diag(table))/sum(table)
也可以使用caret包中的train函数来建模并自动实施10重交叉检验
fitcontrol<-trainControl(method="repeatedcv",number=10,repeats = 3)
tunedf<-data.frame(.cp=seq(0.001,0.1,length.out=10))
treemodel<-train(x=procdata[,-9],y=procdata[,9],method="rpart",trControl=fitcontrol,tuneGrid=tunedf)
treemodel
plot(treemodel)
从上图可以看出，CP参数在0.034可以得到最优的预测准确率
对于目标变量是二元分类变量的情况，除了使用混淆矩阵评价模型以外，还可以使用ROC曲线
很多模型的预测输出允许概率值，混淆矩阵是以0.5为临界点，以划分哪些预测为pos，哪些预测为neg，同时可以相应算出TPR灵敏度和TNR特异度。
除了分类器的训练参数，临界点的选择，也会大大的影响TPR和TNR有时也可以根据具体问题和需要，来选择具体的临界点
如果我们选择一系列的临界点，就会得到一系列的TPR和TNR，将这些值对应的点连接起来，就构成了ROC曲线
ROC曲线可以帮助我们清楚的了解到某个分类器的性能表现，还能方便比较不同分类器的性能
library(ggplot2)
library(caret)
pre<-predict(treemodel,type="prob")
#预测概率prob和实际结果
preObs<-data.frame(prob=pre$pos,obs=procdata$class)
#按预测概率从低到高排序
preObs<-preObs[order(preObs$prob),]
n<-nrow(data)
tpr<-fpr<-rep(0,n)
#根据不同的临界值threshold来计算TPR和FPR，之后绘制成图
for(i in 1:n){
threshold<-preObs$prob[i]
tp<-sum(preObs$prob>threshold&preObs$obs=="pos")
fp<-sum(preObs$prob>threshold&preObs$obs=="neg")
tn<-sum(preObs$prob<threshold&preObs$obs=="neg")
fn<-sum(preObs$prob<threshold&preObs$obs=="pos")
tpr[i]<-tp/(tp+fn) #真正率
fpr[i]<-fp/(tn+fp) #假正率
}
plot(fpr,tpr,type="l")
ROC曲线也可用专门绘制的函数，不仅可以用来画ROC曲线图，还能计算ROC曲线下面积AUC，以评价分类器的综合性嫩，该数值取0-1之间，越大越好，下面的AUC为0.743，要注意这是训练集的结果，不要过于乐观。
library(pROC)
pre<-predict(treemodel,type="prob")$pos
modelroc<-roc(procdata$class,pre)
plot(modelroc,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=TRUE,auc.polygon.col="skyblue",print.thres=TRUE)
还可以使用BP神经网络
library(nnet)
fitcontrol<-trainControl(method="repeatedcv",number=10,repeats = 3)
tunedf<-expand.grid(.decay=0.1,.size=5:10,.bag=T)
nnetmodel<-train(class~.,data=procdata,method="avNNet",trControl=fitcontrol,trace=F,linout=F,tuneGrid=tunedf)
plot(nnetmodel)
网友评论