1、Crime&Shock数据集分析展示,探索性数据分析,Communities and CrimeUnnormalized Data Set （Data source）,2,一、数据预处理,#导入crime数据，修改变量名称，并查看数据属性crim=read.table(crim.txt,sep=,na.string=?)name=read.table(attr_vol.txt)name=name,2colnames(crim)-namesummary(crim) dim(crim),3,观测值：2215 变量数：147部分数据严重缺失,四、犯罪率分布情况,4,在3月份的行业销售旺季，东北地区及
2、北部地区销售额占到公司全月总额的70%，西部地区仅为10%，西部死去市场潜力还需深度挖掘。,可看出violentPerPop nonViolPerPop都出现了不同程度的拖尾特征，考虑对数据进行对数变换,由图可知,四、对数化变换,5,做变换后两变量数据较为对称,由图可知,四、犯罪率地区差异,6,三个地区犯罪率的中位数由西向东递减，分布相对集中，但东部地区出现了较为明显的离群值,7,缺失值处理,nrow(crim!complete.cases(crim),) #缺失值项的总行数#基本每行都有缺失值na.sta=c()for(i in 1:2215)na.sta=c(na.sta,length(w
3、hich(is.na(crimi,)max(na.sta)#缺失值基本在20左右，没有#缺失过于严重的样本，无需删除,8,缺失值处理：临近值插补,从数据集中选取若干条其他属性和它相似的样本(即和它在空间欧式距离最短的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 个样品的空间距离,9,缺失值处理：临近值插补,
4、#先处理非因子型变量的缺失值,需要将以下步骤进行两次for(r in which(!complete.cases(crim1)crim1r,which(is.na(crim1r,) 0.7）,Linear Regression,for(i in 1:5)crim.test=data.frame(crim1ddi,)crim.train=data.frame(crim1-ddi,)lm.vio=lm(violentPerPop.,data=crim.train,c(6:129,146)pre.train=predict(lm.vio)pre.test=predict(lm.vio, crim.t
5、est,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=me
6、an(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),16,Linear Regression,17,Stepwise,#逐步回归：全部、前向、后向lm1both-step(lm.vio,direction= both)lm1forward-step(lm.vio,direction= forward)lm1back-step(lm.
7、vio,direction= backward)final.lm-step(lm.vio)summary(final.lm),18,Stepwise,19,Stepwise,20,Stepwise：diagnose,21,Stepwise,22,Stepwise,23,Conclusion,由结果可以看出，逐步回归调整后的的R-Square为0.6957，模型检验的结果显示，回归的残差项并不满足正态性假定，同时模型的AIC值25573.06依然过高，这启发我们建立更加精确的预测模型。,24,Ridge,25,Ridge,cor(crim,6:129,use=complete.obs) #变量间
8、的相关性，计算时不考虑缺失值symnum(cor(crim,6:129,use=complete.obs) #简单明显示出变量间的相关#只截取部分结果显示，可以看出#变量之间的共线性较为明显,26,Ridge,library(MASS)ridgelm-lm.ridge(violentPerPop.,data=crim.train,c(6:129,146),lambda=seq(0,200,0.01), model =TRUE)names(ridgelm)ridgelm$lambdawhich.min(ridgelm$GCV) #找到GCV 最小时对应的lambda#广义的交叉验证准则GCV，越
9、小越好ridgelm$coef,which.min(ridgelm$GCV) #找到GCV 最小时对应的系数,27,Ridge,#lamda 同GCV 之间关系的图形plot(ridgelm$lambda,ridgelm$GCV,type=l)abline(v=ridgelm$lambdawhich.min(ridgelm$GCV),col=green),28,Lasso,LASSO方法 在线性模型中，人们必须选择合适的变量；比如常用的逐步回归法就是选择显著的变量而抛弃那些不显著的。Tibshirani(1996)1提出了一个新的方法来处理变量选择的问题。该方法在模型系数绝对值的和小于某常数的
10、条件下，谋求残差平方和最小。该方法既提供了如子集选择方法那样的可以解释的模型，也具有岭回归那样的稳定性。它不删除变量，但使得一些回归系数收缩、变小，甚至为0。因而，该方法被称为lasso(least absolute shrinkage and selection operator，最小绝对值收缩和选择算子2)。,29,Lasso,#将数据集中的因子变量gangUnit转换为哑元变量#训练集gangUnit5=rep(0,1772);gangUnit5which(crim.train$gangUnit=5)=1gangUnit10=rep(0,1772);gangUnit10which(cri
11、m.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);gangUnit5which(crim.test$gangUnit=5)=1gangUnit10=rep(0,443);gangUnit10which(crim.test$gangUnit=10)=1crim.test1=cri
12、m.test,c(6:129);crim.test1$gangUnit=NULLcrim.test2=data.frame(crim.test1,1:121,gangUnit5,gangUnit10,crim.test1,-(1:121),30,Lasso(lars),library(lars)attach(crim.train)lassolm-lars(as.matrix(crim.train2),violentPerPop)Lassolm,31,Lasso(msgps),library(msgps)fit - msgps(as.matrix(crim.train2),crim.train$
13、violentPerPop)summary(fit),32,Lasso(lars),lassolm$Cp,33,Lasso(msgps),coef(fit) #extract coefficients at t selected by model selection criteriacbind(coef(fit),ridgelm$coef,which.min(ridgelm$GCV)1:10,)注：如图，根据前三个准则得出的最终模型系数基本一致，最后一列为岭回归Cp原则的系数，可以看出lasso方法可以使更多的系数“收缩”到0。,34,Lasso(msgps),plot(fit,criteri
14、on=cp,main=lasso),35,Lasso(lars),36,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(
15、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 折交叉验证),37,Lasso(lars),pr
16、edict11-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)#
17、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 准则),38,Lasso,# 下面我们来看看上面（根据10 折交叉验证得出的lasso回归方程）它的残差项的正态性检验。shapiro.test(crim.train$
18、violentPerPop-predict11$fit) #训练集上的残差检验#W = 0.8679, p-value 2.2e-16shapiro.test(crim.test$violentPerPop-predict1$fit) #测试集上的残差检验W = 0.8839, p-value 1，最终系数不易收缩到0，这不助于变量的筛选；而对于r1/2时也能得到同样的效果。这一方法提高了原系数估计的稳定性。,42,Generalized elastic net,#elastic netfit2 - msgps(as.matrix(crim.train2),crim.train$violent
19、PerPop,penalty=enet,alpha=0.5)summary(fit2) coef(fit2)plot(fit2,criterion=cp),43,Generalized elastic net,44,Generalized elastic net,#根据Cp原则NMSE.train=c()NMSE.test=c()for(i in 1:5)crim.test=crim1ddi,crim.train=crim1-ddi,#训练集gangUnit5=rep(0,1772);gangUnit5which(crim.train$gangUnit=5)=1gangUnit10=rep(0
20、,1772);gangUnit10which(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);gangUnit5which(crim.test$gangUnit=5)=1gangUnit10=rep(0,443);gangUnit10which(crim.test$ga
21、ngUnit=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),45,Generalized elastic net,fit - msgps(as.matrix(crim.train2),crim.train$violentPerPop,penalty=enet,alpha=0.5)#训练集predict11-predict(fit,as.matrix(crim.tra
22、in2),1NMSE.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),1NMSE.test=c(NMSE.test,mean(crim.test$violentPerPop-as.numeric(predict1)2)/mean(mean(crim.test$violentPerP
23、op)-crim.test$violentPerPop)2)NMSE.train;mean(NMSE.train)NMSE.test;mean(NMSE.test),46,Generalized elastic net,penalty=“genet”(penalty=“alasso”时计算出现奇异值系统报错),47,Model buildingData Mining,48,K临近回归,library(kknn)NMSE=c()MSE=c()for(i in 1:5)crim.test=crim1ddi,crim.train=crim1-ddi,knn1lm-kknn(murders.,crim
24、.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),49,虽然该算法主要用于分类。不用于拟合模型，但其较为稳定，可先利用该模型观察拟合效果以及预测精度再试图从不稳定模型中
25、得到提升。,K临近回归,50,最终达到的预测精度如下，可见该方法的预测精度仍然较低：,回归树,library(rpart) #调用rpart 包crim.test=data.frame(crim1dd1,)crim.train=data.frame(crim1-dd1,)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) #在
26、树中显示分枝的信息。printcp(rt.train) #显示回归树rt.train 每一步得出的sub-trees 的详细信息,51,回归树,52,回归树,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)$
27、violentPerPop)2),n)#tree 训练集上的NMSE 为: 0.3664984 treepre.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 为:
28、 0.4828972,53,回归树：五折交叉验证,for(i in 1:5)crim.test=data.frame(crim1ddi,)crim.train=data.frame(crim1-ddi,)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$violentPe
29、rPop-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.tes
30、t)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),54,回归树,55,Boosting(mboost),for(i in 1:5)crim.test=data.frame(crim1ddi,)crim.train=data.frame(crim1-ddi,)boost.train=blackboost(violentPerPop.,control=boost_control(mstop=50),d
31、ata=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-treep
32、re.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),56
33、,Boosting,57,Bagging,Bagging 是Breiman 提出的与Boosting 相似的技术。Bagging 技术的主要思想是给定一弱学习算法和一训练集1 1 ( , ),., ( , ) n n x y x y 。让该学习算法训练多轮，每轮的训练集由从初始的训练集中随机取出的n 个训练例组成，初始训练例在某轮训练集中可以出现多次或根本不出现。训练之后可得到一个预测函数序列：1,., t h h ，最终的预测函数H 对分类问题采用投票方式，对回归问题采用简单平均方法对新示例进行判别。Bagging 与Boosting 的区别在于Bagging 的训练集的选择是随机的，各轮训
34、练集之间相互独立，而 Boosting 的训练集的选择不是独立的，各轮训练集的选择与前面各轮的学习结果有关；Bagging 的各个预测函数没有权重，而Boosting是有权重的；Bagging 的各个预测函数可以并行生成，而Boosting 的各个预测函数只能顺序生成。对于象神经网络这样极为耗时的学习方法，Bagging 可通过并行训练节省大量时间开销。以下我将通过R 语言中的ipred 包运用Bagging 算法对该数据集进行分析研究。,58,Bagging(ipred),for(i in 1:5)crim.test=data.frame(crim1ddi,)crim.train=data.
35、frame(crim1-ddi,)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$viole
36、ntPerPop)-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,N
37、M.test)MSE.train=c(MSE.train,M.train)MSE.test=c(MSE.test,M.test),59,Bagging,60,随机森林（randomForest）,randomforest.violentPerPop - randomForest(violentPerPop .,data=crim.train,c(6:129,146),ntree=500,importance=TRUE,proximity=TRUE)randomforest.violentPerPop$importance #查看解释变量对模型的贡献性的大小randomforest.violen
38、tPerPop$importanceSD,61,随机森林（randomForest）,#贡献度最大的前十个变量names(crim.train,c(6:129,146)sort(randomforest.violentPerPop$importance,1,decreasing=T,index.return=T)$ix1:10plot(randomforest.violentPerPop$importanceSD)identify(1:124,randomforest.violentPerPop$importanceSD,labels=names(randomforest.violentPer
39、Pop$importanceSD),62,随机森林（randomForest）,63,随机森林（randomForest）,64,总结,65,Shock Data(Data source),66,数据预处理,67,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)
40、,数据描述,68,69,数据描述,70,数据描述,71,数据描述,Record不同时的shock type完全相同,72,数据描述,由图可知,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
41、, 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)else if (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),
