R数据挖掘实例

Crime&Shock数据集分析展示探索性数据分析CommunitiesandCrimeUnnormalizedDataSet（Datasource）一、数据预处理#导入crime数据，修改变量名称，并查看数据属性crim=read.table("crim.txt",sep=",",na.string="?")name=read.table("attr_vol.txt")name=name[,2]colnames(crim)<-namesummary(crim)dim(crim)观测值：2215变量数：147部分数据严重缺失四、犯罪率分布情况在3月份的行业销售旺季，东北地区及北部地区销售额占到公司全月总额的70%，西部地区仅为10%，西部死去市场潜力还需深度挖掘。可看出violentPerPopnonViolPerPop都出现了不同程度的拖尾特征，考虑对数据进行对数变换由图可知四、对数化变换做变换后两变量数据较为对称由图可知四、犯罪率地区差异三个地区犯罪率的中位数由西向东递减，分布相对集中，但东部地区出现了较为明显的离群值缺失值处理nrow(crim[!complete.cases(crim),])##缺失值项的总行数#基本每行都有缺失值na.sta=c()for(iin1:2215){na.sta=c(na.sta,length(which(is.na(crim[i,]))))}max(na.sta)#缺失值基本在20左右，没有#缺失过于严重的样本，无需删除缺失值处理：临近值插补从数据集中选取若干条其他属性和它相似的样本(即和它在空间欧式距离最短的n条样本)，求其中位数进行插补crim$gangUnit=as.factor(gangUnit)crim1=crim[,c(6:126,128:147)]library(cluster)##调用R语言中的cluster包#计算空间距离dist.mtx<-as.matrix(daisy(crim1,stand=T))##计算这2215个样品的空间距离缺失值处理：临近值插补#先处理非因子型变量的缺失值,需要将以下步骤进行两次for(rinwhich(!complete.cases(crim1)))crim1[r,which(is.na(crim1[r,]))]<-apply(data.frame(crim1[c(as.integer(names(sort(dist.mtx[r,])[2:20]))),which(is.na(crim1[r,]))]),2,median,na.rm=T)#再处理因子型变量的for(rin1:2215){if(is.na(gangUnit[r])){index=sort(dist.mtx[r,],index.return=T)$ixif(all(is.na(gangUnit[index[2:20]])))gangUnit[r]=gangUnit[intersect(index,which(is.na(gangUnit)==F))][1]else{gangUnit[r]=levels(gangUnit[index[2:11]])[which.max(table(gangUnit[index[2:11]]))]}}}crim2=data.frame(cbind(crim[,1:5],crim1[,1:126],gangUnit,crim1[,-(1:126)]))Crossvalidation#设置五折交叉验证n=2215;zz1=1:n#zz1为所有观测值(行)的下标zz2=rep(1:5,ceiling(n/5))[1:n]set.seed(200);zz2=sample(zz2,n)#zz2为1:5的随机排列zz2[1:100]#dd保存每一折下标，令testset选其中之一，则共可做五次交叉验证dd=list()for(iin1:5)dd[[i]]=zz1[zz2==i]ddModelbuilding——traditionalmethods LinearRegressioncrim.test=crim1[dd[[1]],]crim.train=crim1[-dd[[1]],]lm.vio=lm(violentPerPop~.,data=crim.train[,c(6:129,146)])summary(lm.vio)#对因变量做对数变换lm.logvio=lm(log(violentPerPop)~.,data=crim.train[,c(6:129,146)])summary(lm.logvio)shapiro.test(lm.vio$residuals)shapiro.test(lm.logvio$residuals)LinearRegression结果显示对因变量做对数变换后并未使模型显著度明显增强。注：由于空间原因只显示未对数化处理模型的部分结果LinearRegression：anova注：由于空间原因只显示部分结果LinearRegression：anovaav=anova(lm.vio)#显示回归中系数不显著的前十个变量名称names(crim[,6:129])[sort(av[,5],decreasing=T,index.return=T)$ix[1:20]]注：红线标注的为两个模型中都非常不显著变量（pvalue>0.7）LinearRegressionfor(iin1:5){crim.test=data.frame(crim1[dd[[i]],])crim.train=data.frame(crim1[-dd[[i]],])lm.vio=lm(violentPerPop~.,data=crim.train[,c(6:129,146)])pre.train=predict(lm.vio)pre.test=predict(lm.vio,crim.test[,c(6:129,146)])NM.train=mean((crim.train$violentPerPop-pre.train)^2)/mean((mean(crim.train$violentPerPop)-crim.train$violentPerPop)^2)NM.test=mean((crim.test$violentPerPop-pre.test)^2)/mean((mean(crim.test$violentPerPop)-crim.test$violentPerPop)^2)M.train=mean((crim.train$violentPerPop-pre.train)^2)M.test=mean((crim.test$violentPerPop-pre.test)^2)NMSE.train=c(NMSE.train,NM.train)NMSE.test=c(NMSE.test,NM.test)MSE.train=c(MSE.train,M.train)MSE.test=c(MSE.test,M.test)}LinearRegressionStepwise#逐步回归：全部、前向、后向lm1both<-step(lm.vio,direction="both")lm1forward<-step(lm.vio,direction="forward")lm1back<-step(lm.vio,direction="backward")final.lm<-step(lm.vio)summary(final.lm)StepwiseStepwiseStepwise：diagnoseStepwiseStepwiseConclusion由结果可以看出，逐步回归调整后的的R-Square为0.6957，模型检验的结果显示，回归的残差项并不满足正态性假定，同时模型的AIC值25573.06依然过高，这启发我们建立更加精确的预测模型。RidgeRidgecor(crim[,6:129],use="complete.obs")#变量间的相关性，计算时不考虑缺失值symnum(cor(crim[,6:129],use="complete.obs"))#简单明显示出变量间的相关#只截取部分结果显示，可以看出#变量之间的共线性较为明显Ridgelibrary(MASS)ridgelm<-lm.ridge(violentPerPop~.,data=crim.train[,c(6:129,146)],lambda=seq(0,200,0.01),model=TRUE)names(ridgelm)ridgelm$lambda[which.min(ridgelm$GCV)]##找到GCV最小时对应的lambda##广义的交叉验证准则—GCV，越小越好ridgelm$coef[,which.min(ridgelm$GCV)]#找到GCV最小时对应的系数Ridge#lamda同GCV之间关系的图形plot(ridgelm$lambda,ridgelm$GCV,type="l")abline(v=ridgelm$lambda[which.min(ridgelm$GCV)],col="green")LassoLASSO方法在线性模型中，人们必须选择合适的变量；比如常用的逐步回归法就是选择显著的变量而抛弃那些不显著的。Tibshirani(1996)[1]提出了一个新的方法来处理变量选择的问题。该方法在模型系数绝对值的和小于某常数的条件下，谋求残差平方和最小。该方法既提供了如子集选择方法那样的可以解释的模型，也具有岭回归那样的稳定性。它不删除变量，但使得一些回归系数收缩、变小，甚至为0。因而，该方法被称为lasso(leastabsoluteshrinkageandselectionoperator，最小绝对值收缩和选择算子[2])。Lasso#将数据集中的因子变量gangUnit转换为哑元变量#训练集gangUnit5=rep(0,1772);gangUnit5[which(crim.train$gangUnit==5)]=1gangUnit10=rep(0,1772);gangUnit10[which(crim.train$gangUnit==10)]=1crim.train1=crim.train[,c(6:129)];crim.train1$gangUnit=NULLcrim.train2=data.frame(crim.train1[,1:121],gangUnit5,gangUnit10,crim.train1[,-(1:121)])#测试集gangUnit5=rep(0,443);gangUnit5[which(crim.test$gangUnit==5)]=1gangUnit10=rep(0,443);gangUnit10[which(crim.test$gangUnit==10)]=1crim.test1=crim.test[,c(6:129)];crim.test1$gangUnit=NULLcrim.test2=data.frame(crim.test1[,1:121],gangUnit5,gangUnit10,crim.test1[,-(1:121)])Lasso(lars)library(lars)attach(crim.train)lassolm<-lars(as.matrix(crim.train2),violentPerPop)LassolmLasso(msgps)library(msgps)fit<-msgps(as.matrix(crim.train2),crim.train$violentPerPop)summary(fit)Lasso(lars)lassolm$CpLasso(msgps)coef(fit)#extractcoefficientsattselectedbymodelselectioncriteriacbind(coef(fit),ridgelm$coef[,which.min(ridgelm$GCV)])[1:10,])注：如图，根据前三个准则得出的最终模型系数基本一致，最后一列为岭回归Cp原则的系数，可以看出lasso方法可以使更多的系数“收缩”到0。Lasso(msgps)plot(fit,criterion="cp",main="lasso")Lasso(lars)Lasso(lars)predict11<-predict(lassolm,data.matrix(crim.train2),mode="fraction",s=best,type="fit")predict1<-predict(lassolm,data.matrix(crim.test2),mode="fraction",s=best,type="fit")cat("lasso训练集上的NMSE为:",mean((crim.train$violentPerPop-as.numeric(predict11$fit))^2)/mean((mean(crim.train$violentPerPop)-crim.train$violentPerPop)^2),"\n")lasso训练集上的NMSE为:0.3050613(依据—10折交叉验证)cat("lasso测试集上的NMSE为:",mean((crim.test$violentPerPop-as.numeric(predict1$fit))^2)/mean((mean(crim.test$violentPerPop)-crim.test$violentPerPop)^2),"\n")lasso测试集上的NMSE为:0.3388574(依据—10折交叉验证)Lasso(lars)predict11<-predict(lassolm,data.matrix(crim.train2),mode="step",s=best,type="fit")predict1<-predict(lassolm,data.matrix(crim.test2),mode="step",s=best,type="fit")cat("lasso训练集上的NMSE为:",mean((crim.train$violentPerPop-as.numeric(predict11$fit))^2)/mean((mean(crim.train$violentPerPop)-crim.train$violentPerPop)^2),"\n")#lasso训练集上的NMSE为:0.3248435(依据—Cp准则)cat("lasso测试集上的NMSE为:",mean((crim.test$violentPerPop-as.numeric(predict1$fit))^2)/mean((mean(crim.test$violentPerPop)-crim.test$violentPerPop)^2),"\n")#lasso测试集上的NMSE为:0.3352093(依据—Cp准则)Lasso##下面我们来看看上面（根据10折交叉验证得出的lasso回归方程）它的残差项的正态性检验。shapiro.test(crim.train$violentPerPop-predict11$fit)##训练集上的残差检验#W=0.8679,p-value<2.2e-16shapiro.test(crim.test$violentPerPop-predict1$fit)##测试集上的残差检验W=0.8839,p-value<2.2e-16可以看出上述模型的残差项仍然不满足正态性的分布假定，虽然达到了较好的预测精度，但是我们的模型还需要进一步改进Ridge&Lasso实际上，岭回归和lasso方法都是用于处理模型变量共线性的方式。在传统线性回归当中，变量共线性极有可能导致估计的回归系数的不稳定（这一特性在样本总量较少而变量总数较多的情形下尤其显著），这一缺点可以通过对回归系数加以约束得到改善。不妨设约束函数为如下形式，则岭回归和lasso回归可以在参数族的形式下得到统一（岭回归：r=2;lasso:r=1）虽然只是参数取法上的不一致，但由于约束函数凹凸函数不同的特性（r<1为凹函数，r>1为凸函数）导致了最终系数特点的不一致。对于岭回归，随着t的减小，模型系数绝对值一般“步调一致”的收缩；对于lasso回归，则更倾向于将部分变量的系数收敛到0而保持另外一些变量系数的大小来达到约束条件（这实际上增加了最终系数的差异）对于r取值的其他情形而导致的后果也可由函数凹凸性作相似推断。具体参见FastSparseRegressionandClassificationJeromeH.Friedman2008Ridge&Lasso这一参数族中r的选取也带来很多问题，对于r>1，最终系数不易收缩到0，这不助于变量的筛选；而对于r<1，随着t的减小（可以一一对应到拉格朗日方法中的λ）其系数变化呈现不连续的特征，这一问题可以通过加入“弹性系数”的方式（Generalizedelasticnet）得以改进（在一定的参数取值条件下）。Generalizedelasticnetmsgps这一软件包给我们提供了解决 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 。根据Friedman论文中结论，在enet方法中，当α在(0,1)中时，其得到的系数变化较为连续（α=0对应lasso方法）；而在genet方法中α>1/2时也能得到同样的效果。这一方法提高了原系数估计的稳定性。Generalizedelasticnet#elasticnetfit2<-msgps(as.matrix(crim.train2),crim.train$violentPerPop,penalty="enet",alpha=0.5)summary(fit2)coef(fit2)plot(fit2,criterion="cp")GeneralizedelasticnetGeneralizedelasticnet#根据Cp原则NMSE.train=c()NMSE.test=c()for(iin1:5){crim.test=crim1[dd[[i]],]crim.train=crim1[-dd[[i]],]#训练集gangUnit5=rep(0,1772);gangUnit5[which(crim.train$gangUnit==5)]=1gangUnit10=rep(0,1772);gangUnit10[which(crim.train$gangUnit==10)]=1crim.train1=crim.train[,c(6:129)];crim.train1$gangUnit=NULLcrim.train2=data.frame(crim.train1[,1:121],gangUnit5,gangUnit10,crim.train1[,-(1:121)])#测试集gangUnit5=rep(0,443);gangUnit5[which(crim.test$gangUnit==5)]=1gangUnit10=rep(0,443);gangUnit10[which(crim.test$gangUnit==10)]=1crim.test1=crim.test[,c(6:129)];crim.test1$gangUnit=NULLcrim.test2=data.frame(crim.test1[,1:121],gangUnit5,gangUnit10,crim.test1[,-(1:121)])Generalizedelasticnetfit<-msgps(as.matrix(crim.train2),crim.train$violentPerPop,penalty="enet",alpha=0.5)#训练集predict11<-predict(fit,as.matrix(crim.train2))[,1]NMSE.train=c(NMSE.train,mean((crim.train$violentPerPop-as.numeric(predict11))^2)/mean((mean(crim.train$violentPerPop)-crim.train$violentPerPop)^2))#测试集predict1<-predict(fit,as.matrix(crim.test2))[,1]NMSE.test=c(NMSE.test,mean((crim.test$violentPerPop-as.numeric(predict1))^2)/mean((mean(crim.test$violentPerPop)-crim.test$violentPerPop)^2))}NMSE.train;mean(NMSE.train)NMSE.test;mean(NMSE.test)Generalizedelasticnetpenalty=“genet”(penalty=“alasso”时计算出现奇异值系统报错)Modelbuilding——DataMiningK临近回归library(kknn)NMSE=c()MSE=c()for(iin1:5){crim.test=crim1[dd[[i]],]crim.train=crim1[-dd[[i]],]knn1lm<-kknn(murders~.,crim.train[,6:130],crim.test[,6:130],k=1,distance=1,kernel="rectangular")NM=mean((crim.test$murders-knn1lm$fitted.values)^2)/mean((mean(crim.test$murders)-crim.test$murders)^2)M=mean((crim.test$murders-knn1lm$fitted.values)^2)NMSE=c(NMSE,NM)MSE=c(MSE,M)}虽然该算法主要用于分类。不用于拟合模型，但其较为稳定，可先利用该模型观察拟合效果以及预测精度再试图从不稳定模型中得到提升。K临近回归最终达到的预测精度如下，可见该方法的预测精度仍然较低：回归树library(rpart)##调用rpart包crim.test=data.frame(crim1[dd[[1]],])crim.train=data.frame(crim1[-dd[[1]],])rt.train=rpart(violentPerPop~.,data=crim.train[,c(6:129,146)])plot(rt.train,uniform=T,branch=1,margin=0.1,cex=0.9,main="violentPerPop")##画出回归树text(rt.train,cex=0.75)##在树中显示分枝的信息。printcp(rt.train)##显示回归树rt.train每一步得出的sub-trees的详细信息回归树回归树treepre.train=predict(rt.train,crim.train[,c(6:129,146)])cat("tree训练集上的NMSE为:",mean((crim.train[,c(6:129,146)]$violentPerPop-as.numeric(treepre.train))^2)/mean((mean(crim.train[,c(6:129,146)]$violentPerPop)-crim.train[,c(6:129,146)]$violentPerPop)^2),"\n")#tree训练集上的NMSE为:0.3664984treepre.test=predict(rt.train,crim.test[,c(6:129,146)])cat("tree训练集上的NMSE为:",mean((crim.test[,c(6:129,146)]$violentPerPop-as.numeric(treepre.test))^2)/mean((mean(crim.test[,c(6:129,146)]$violentPerPop)-crim.test[,c(6:129,146)]$violentPerPop)^2),"\n")#tree测试集上的NMSE为:0.4828972回归树：五折交叉验证for(iin1:5){crim.test=data.frame(crim1[dd[[i]],])crim.train=data.frame(crim1[-dd[[i]],])rt.train=rpart(violentPerPop~.,data=crim.train[,c(6:129,146)])treepre.train=predict(rt.train,crim.train[,c(6:129,146)])treepre.test=predict(rt.train,crim.test[,c(6:129,146)])NM.train=mean((crim.train$violentPerPop-treepre.train)^2)/mean((mean(crim.train$violentPerPop)-crim.train$violentPerPop)^2)NM.test=mean((crim.test$violentPerPop-treepre.test)^2)/mean((mean(crim.test$violentPerPop)-crim.test$violentPerPop)^2)M.train=mean((crim.train$violentPerPop-treepre.train)^2)M.test=mean((crim.test$violentPerPop-treepre.test)^2)NMSE.train=c(NMSE.train,NM.train)NMSE.test=c(NMSE.test,NM.test)MSE.train=c(MSE.train,M.train)MSE.test=c(MSE.test,M.test)}回归树Boosting(mboost)for(iin1:5){crim.test=data.frame(crim1[dd[[i]],])crim.train=data.frame(crim1[-dd[[i]],])boost.train=blackboost(violentPerPop~.,control=boost_control(mstop=50),data=crim.train[,c(6:129,146)])treepre.train=predict(boost.train,crim.train[,c(6:129,146)])treepre.test=predict(boost.train,crim.test[,c(6:129,146)])NM.train=mean((crim.train$violentPerPop-treepre.train)^2)/mean((mean(crim.train$violentPerPop)-crim.train$violentPerPop)^2)NM.test=mean((crim.test$violentPerPop-treepre.test)^2)/mean((mean(crim.test$violentPerPop)-crim.test$violentPerPop)^2)M.train=mean((crim.train$violentPerPop-treepre.train)^2)M.test=mean((crim.test$violentPerPop-treepre.test)^2)NMSE.train=c(NMSE.train,NM.train)NMSE.test=c(NMSE.test,NM.test)MSE.train=c(MSE.train,M.train)MSE.test=c(MSE.test,M.test)}}BoostingBaggingBagging是Breiman提出的与Boosting相似的技术。Bagging技术的主要 思想 教师资格思想品德鉴定表下载浅论红楼梦的主题思想员工思想动态调查问卷论语教育思想学生思想教育讲话稿 是给定一弱学习算法和一训练集11(,),...,(,)nnxyxy。让该学习算法训练多轮，每轮的训练集由从初始的训练集中随机取出的n个训练例组成，初始训练例在某轮训练集中可以出现多次或根本不出现。训练之后可得到一个预测函数序列：1,...,thh，最终的预测函数H对分类问题采用投票方式，对回归问题采用简单平均方法对新示例进行判别。Bagging与Boosting的区别在于Bagging的训练集的选择是随机的，各轮训练集之间相互独立，而Boosting的训练集的选择不是独立的，各轮训练集的选择与前面各轮的学习结果有关；Bagging的各个预测函数没有权重，而Boosting是有权重的；Bagging的各个预测函数可以并行生成，而Boosting的各个预测函数只能顺序生成。对于象神经网络这样极为耗时的学习方法，Bagging可通过并行训练节省大量时间开销。以下我将通过R语言中的ipred包运用Bagging算法对该数据集进行分析研究。Bagging(ipred)for(iin1:5){crim.test=data.frame(crim1[dd[[i]],])crim.train=data.frame(crim1[-dd[[i]],])bagging.vio=bagging(violentPerPop~.,data=crim.train[,c(6:129,146)],coob=T,control=rpart.control(xval=10))pre.train=predict(bagging.vio,crim.train[,c(6:129,146)])pre.test=predict(bagging.vio,crim.test[,c(6:129,146)])NM.train=mean((crim.train$violentPerPop-pre.train)^2)/mean((mean(crim.train$violentPerPop)-crim.train$violentPerPop)^2)NM.test=mean((crim.test$violentPerPop-pre.test)^2)/mean((mean(crim.test$violentPerPop)-crim.test$violentPerPop)^2)M.train=mean((crim.train$violentPerPop-pre.train)^2)M.test=mean((crim.test$violentPerPop-pre.test)^2)NMSE.train=c(NMSE.train,NM.train)NMSE.test=c(NMSE.test,NM.test)MSE.train=c(MSE.train,M.train)MSE.test=c(MSE.test,M.test)}Bagging随机森林（randomForest）randomforest.violentPerPop<-randomForest(violentPerPop~.,data=crim.train[,c(6:129,146)],ntree=500,importance=TRUE,proximity=TRUE)randomforest.violentPerPop$importance##查看解释变量对模型的贡献性的大小randomforest.violentPerPop$importanceSD随机森林（randomForest）#贡献度最大的前十个变量names(crim.train[,c(6:129,146)])[sort(randomforest.violentPerPop$importance[,1],decreasing=T,index.return=T)$ix[1:10]]plot(randomforest.violentPerPop$importanceSD)identify(1:124,randomforest.violentPerPop$importanceSD,labels=names(randomforest.violentPerPop$importanceSD))随机森林（randomForest）随机森林（randomForest）总结方法NMSE训练集测试集传统统计方法线性回归0.29000.3729逐步回归0.29640.3774lasso0.30510.3389Generalizedelasticnet0.29380.3677数据挖掘K临近回归无0.6257回归树0.35290.4627Boosting0.27530.3496Bagging0.26950.3589随机森林0.32580.3228ShockData(Datasource)数据预处理shock=read.table("shock.txt",header=T)head(shock)shock$SHOCK_TYP=as.factor(shock$SHOCK_TYP)shock$SURVIVE=as.factor(shock$SURVIVE)shock$SEX=as.factor(shock$SEX)shock$RECORD=as.factor(shock$RECORD)数据描述数据描述数据描述数据描述Record不同时的shocktype完全相同数据描述由图可知Record不同时各观测值有差异但不明显，由箱盒图可以看到HT与RCI中存在离群值贝叶斯分类distinguish.bayes<-function(TrnX,TrnG,p=rep(1,length(levels(TrnG))),TstX=NULL,var.equal=FALSE){if(is.factor(TrnG)==FALSE){mx<-nrow(TrnX);mg<-nrow(TrnG)TrnX<-rbind(TrnX,TrnG)TrnG<-factor(rep(1:2,c(mx,mg)))}if(is.null(TstX)==TRUE)TstX<-TrnXif(is.vector(TstX)==TRUE)TstX<-t(as.matrix(TstX))elseif(is.matrix(TstX)!=TRUE)TstX<-as.matrix(TstX)if(is.matrix(TrnX)!=TRUE)TrnX<-as.matrix(TrnX)nx<-nrow(TstX)blong<-matrix(rep(0,nx),nrow=1,dimnames=list("blong",1:nx))g<-length(levels(TrnG))mu<-matrix(0,nrow=g,ncol=ncol(TrnX))贝叶斯分类for(iin1:g)mu[i,]<-colMeans(TrnX[TrnG==i,])D<-matrix(0,nrow=g,ncol=nx)if(var.equal==TRUE||var.equal==T){for(iin1:g){d2<-mahalanobis(TstX,mu[i,],var(TrnX))D[i,]<-d2-2*log(p[i])}}else{for(iin1:g){S<-var(TrnX[TrnG==i,])d2<-mahalanobis(TstX,mu[i,],S)D[i,]<-d2-2*log(p[i])-log(det(S))}}for(jin1:nx){dmin<-Inffor(iin1:g)if(D[i,j]<dmin){dmin<-D[i,j];blong[j]<-i}}blong}贝叶斯分类shock2=shockshock2$SURVIVE=as.numeric(shock2$SURVIVE)shock2$SEX=as.numeric(shock2$SEX)shock2$RECORD=as.numeric(shock2$RECORD)distinguish.bayes(shock2[,-c(1,6)],shock2$SHOCK_TYP)#将因子数据转化为01变量shock2=shock[,-c(which(names(shock)=="SHOCK_TYP"),which(names(shock)=="SURVIVE"),which(names(shock)=="SEX"),which(names(shock)=="RECORD"))]shock3=cbind(shock2,sex=as.numeric(sex[,1]),survive=as.numeric(survive[,1]),record=as.numeric(record[,1]))决策树library(rpart)##调用rpart包rt.type<-rpart(SHOCK_TYP~.,data=shock[-dd[[1]],-1])rt.type决策树决策树printcp(rt.type)##显示回归树rt.type每一步得出的sub-trees的详细信息决策树#五折交叉验证train=c();test=c();for(iin1:5){rt.type<-rpart(SHOCK_TYP~.,data=shock[-dd[[i]],])c1=t(table(predict(rt.type,shock[-dd[[i]],],type="class"),shock[-dd[[i]],6]))c2=t(table(predict(rt.type,shock[dd[[i]],],type="class"),shock[dd[[i]],6]))train=c(train,sum(diag(c1))/sum(c1))test=c(test,sum(diag(c2))/sum(c2))}train;mean(train)test;mean(test)最邻近算法#五折交叉验证,对k=1:15做循环final.test=c();final.mean=c()for(jin1:15){test=c();for(iin1:5){knn1lm<-kknn(SHOCK_TYP~.,shock[-dd[[i]],],shock[dd[[i]],],k=j,distance=1,kernel="rectangular")c=table(shock[dd[[i]],6],knn1lm$fitted.values)test=c(test,sum(diag(c))/sum(c))}final.test=cbind(final.test,test)final.mean=c(final.mean,mean(test))}final.mean;max(final.mean)最邻近算法boostinglibrary(adabag)a=boosting(SHOCK_TYP~.,data=shock[-dd[[1]],],mfinal=15)aa.pred<-predict.boosting(a,newdata=shock[-dd[[1]],]);a.predboostingbarplot(a$importance)boostingtrain=c();test=c();importance=c();for(iin1:5){a=boosting(SHOCK_TYP~.,data=shock[-dd[[i]],],mfinal=15)a.pred<-predict.boosting(a,newdata=shock[-dd[[i]],])a.predt<-predict.boosting(a,newdata=shock[dd[[i]],])train=c(train,1-a.pred$error)test=c(test,1-a.predt$error)importance=cbind(importance,names(a$importance)[sort(a$importance,decreasing=T,index.return=T)$ix[1:5]])}train;mean(train)test;mean(test)importanceboostingbagginglibrary(adabag)train=c();test=c();importance=c();for(iin1:5){a=bagging(SHOCK_TYP~.,data=shock[-dd[[i]],],mfinal=15)a.pred<-predict.bagging(a,newdata=shock[-dd[[i]],])a.predt<-predict.bagging(a,newdata=shock[dd[[i]],])train=c(train,1-a.pred$error)test=c(test,1-a.predt$error)importance=cbind(importance,names(a$importance)[sort(a$importance,decreasing=T,index.return=T)$ix[1:5]])}train;mean(train)test;mean(test)importancebagging随机森林library(randomForest)w=randomForest(SHOCK_TYP~.,data=shock[-dd[[1]],],importance=TRUE,proximity=TRUE)barplot(w$importance)par(mfrow=c(2,1))barplot(t(importance(w))[7,],cex.names=0.6,main="MeanDecreaseAccuracy")barplot(t(importance(w))[8,],cex.names=0.6,main="MeanDecreaseGini")随机森林随机森林print(w)随机森林#五折交叉验证train=c();test=c();for(iin1:5){w=randomForest(SHOCK_TYP~.,data=shock[-dd[[i]],],importance=TRUE,proximity=TRUE)c1=t(table(predict(w),shock[-dd[[i]],6]))c2=t(table(predict(w,shock[dd[[i]],]),shock[dd[[i]],6]))train=c(train,sum(diag(c1))/sum(c1))test=c(test,sum(diag(c2))/sum(c2))}train;mean(train)test;mean(test)随机森林神经网络（nnet）#实验法找最优节点，最后选取隐结点数为23左右为最优设计train=c();test=c();for(iin5:30){a=nnet(shock[-dd[[1]],-6],b[-dd[[1]],],size=i,rang=0.1,decay=5e-4,maxit=200);c1=test.cl(b[-dd[[1]],],predict(a,shock[-dd[[1]],-6]))c2=test.cl(b[dd[[1]],],predict(a,shock[dd[[1]],-6]))train=c(train,sum(diag(c1))/sum(c1))test=c(test,sum(diag(c2))/sum(c2))}result=data.frame(size=5:30,train,test)神经网络（nnet）#五折交叉验证train=c();test=c();for(iin1:5){a=nnet(shock[-dd[[i]],-6],b[-dd[[i]],],size=23,rang=0.1,decay=5e-4,maxit=200);c1=test.cl(b[-dd[[i]],],predict(a,shock[-dd[[i]],-6]))c2=test.cl(b[dd[[i]],],predict(a,shock[dd[[i]],-6]))train=c(train,sum(diag(c1))/sum(c1))test=c(test,sum(diag(c2))/sum(c2))}train;mean(train)test;mean(test)神经网络（nnet）注：左图为size=23、25、27的实验情形，推断size=27出现过拟合现象，同时实验结果较为不稳定支持向量机#五折交叉验证library(e1071)train=c();test=c()for(iin1:5){model<-svm(SHOCK_TYP~.,data=shock[-dd[[i]],],kernal="sigmoid")pred.train<-fitted(model)r1=table(pred.train,shock$SHOCK_TYP[-dd[[i]]])pred.test<-predict(model,shock[dd[[i]],-6])r2=table(pred.test,shock$SHOCK_TYP[dd[[i]]])train=c(train,sum(diag(r1))/sum(r1))test=c(test,sum(diag(r2))/sum(r2))}train;mean(train)test;mean(test)支持向量机（e1071）总结方法准确率训练集测试集决策树0.60840.3360最邻近算法无0.4424Boosting0.72350.4240Bagging0.71020.4154随机森林0.53210.5441神经网络(size=25)0.44140.3275支持向量机0.66820.4335Thanks