现代回归和分类

现代回归和分类:

算法建模的回归和分类吴喜之算法建模的回归和分类经典的回归和分类（判别）模型是可以写成公式的。但是另外一些回归和分类的方法是体现在算法之中，其具体形式是计算机程序。广义地说，算法模型实际上包含了经典模型如果说，起源于前计算机时代的经典统计目前大大受惠于计算机的发展，那么，没有计算机，就不可能存在算法建模。在处理巨大的数据集上，在对付称为维数诅咒的大量的变量上，在不用假定总体分布的情况时，在对付众多竞争模型方面，算法建模较经典建模有很多不可比拟的优越性。决策树：分类树和回归树

例(数据shuttle.txt)

library(MASS);shuttle[1:10,]这个数据是关于美国航天飞机在各种条件下是否自动着陆的决策问题[1]。有256行及7列。头六行为作为自变量的定性变量，而最后一列为因变量。自变量包括稳定性(stability，取值stab/xstab)、误差大小(error，取值(MM/SS/LX/XL)、信号(sign，取值pp/nn)、风向(wind，取值head/tail)、风力(magn，取值(Light/Medium/Strong/Out)、能见度(vis，取值yes/no)，因变量为是否用自动着陆系统(use，取值auto/noauto)。

[1]

数据源于D.Michie(1989)Problemsofcomputer-aidedconceptformation.InApplicationsofExpertSystems2,ed.J.R.Quinlan,TuringInstitutePress/Addison-Wesley,pp.310–333.例(数据shuttle.txt).library(MASS);shuttle[1:10,]m=256;set.seed(2);samp=sample(1:m,floor(m/10));tsamp=setdiff(1:m,samp)library(rpart.plot);(b=rpart(use~.,shuttle,subset=tsamp));b;plot(b);text(b,use.n=T)t(table(predict(b,shuttle[tsamp,],type="class"),shuttle[tsamp,7]))rpart.plot(b,type=4,extra=6,faclen=T);rpart.plot(b,type=0,extra=6,faclen=T);rpart.plot(b,type=1,extra=6,faclen=T);rpart.plot(b,type=2,extra=6,faclen=T);例10.1(数据shuttle.txt).t(table(predict(b,shuttle[tsamp,],type="class"),shuttle[tsamp,7]))t(table(predict(b,shuttle[samp,],type="class"),shuttle[samp,7]))kyphosisdataThekyphosisdataframehas81rowsand4columns.representingdataonchildrenwhohavehadcorrectivespinalsurgery

Kyphosis:afactorwithlevelsabsentpresentindicatingifakyphosis(atypeofdeformation)waspresentaftertheoperation.

Age:inmonths

Number:thenumberofvertebraeinvolved

Start:thenumberofthefirst(topmost)vertebraoperatedon.

library(rpart.plot)fit<-rpart(Kyphosis~Age+Number+Start,data=kyphosis)fit2<-rpart(Kyphosis~Age+Number+Start,data=kyphosis,parms=list(prior=c(.65,.35),split='information'))fit3<-rpart(Kyphosis~Age+Number+Start,data=kyphosis,control=rpart.control(cp=.05))par(mfrow=c(1,3),xpd=NA);plot(fit);text(fit,use.n=TRUE)plot(fit2);text(fit2,use.n=TRUE);plot(fit3);text(fit3,use.n=TRUE)剪枝和画图par(mfrow=c(1,3),xpd=NA);rpart.plot(fit,type=2,extra=6)rpart.plot(fit2,type=2,extra=6);rpart.plot(fit3,type=2,extra=6);par(mfrow=c(1,1))library(rpart)fit<-rpart(Kyphosis~Age+Number+Start,data=kyphosis)predict(fit,type="prob")#classprobabilities(default)predict(fit,type="vector")#levelnumberspredict(fit,type="class")#factorpredict(fit,type="matrix")#levelnumber,classfrequencies,probabilities预测library(rpart)library(rpart.plot)data(kyphosis)kyphosis.rp<-rpart(Kyphosis~Age+Number+Start,data=kyphosis,subset=1:70)kyphosis.rp;plot(kyphosis.rp);text(kyphosis.rp,use.n=T)预测(2)rpart.plot(kyphosis.rp,type=2,extra=6)library(rpart)kyphosis1<-kyphosis[71:81,]predict(kyphosis.rp,kyphosis1,type="class")table(predict(kyphosis.rp,kyphosis1,type="class"),kyphosis[71:81,1])

预测(2)例10.2(例9.5数据iris.txt).

library(MASS);m=150;set.seed(10)samp<-c(sample(1:50,25),sample(51:100,25),sample(101:150,25));tsamp=setdiff(1:m,samp);library(rpart.plot)(b=rpart(Species~.,iris,subset=tsamp));plot(b);text(b,use.n=T)rpart.plot(b,type=2,extra=6)IRISt(table(predict(b,iris[tsamp,],type="class"),iris[tsamp,5]))t(table(predict(b,iris[samp,],type="class"),iris[samp,5]))(数据wine.txt)这是关于意大利一个地区的葡萄酒数据[1]，该数据是对该地区的三种不同培育品种的葡萄所酿造的酒的13种要素的化学分析结果，一共有178个观测值。我们希望用这些数据来建立一个模型，利用这些要素的特性来判断是那个品种的葡萄所酿造的酒。因变量是Class(品种)，有1、2、3个哑元取值，而13个自变量为：Alcohol(酒精)、Malic.acid(苹果酸)、Ash(灰分)、Alcalinityofash(灰的碱性)、Magnesium(镁)、Totalphenols(苯酚总量)、Flavanoids[2]、Nonflavanoidphenols(非flavanoid苯酚)、Proanthocyanins[3]、Colorintensity(颜色强度)、Hue(色调)、OD280/OD315ofdilutedwines(稀释的酒的蛋白质浓度的光谱度量)、Proline(脯氨酸)。根据这13个自变量，我们建立了下面的一个决策树：

w=read.table("f:/adbook/data/wine1.txt",header=T,sep=",")w$Class=factor(w$Class)library(rpart)(fit<-rpart(Class~Flavanoids+Censity+Proline,data=w))plot(fit);text(fit)rpart.plot(fit,type=2,extra=6)思考一下：无论自变量是连续的还是定性变量，分类树的原理都一样。连续变量在根据需要离散化之后和分类变量就类似了。后面要介绍的回归树的因变量也离散化了（但也不是事先离散化），于是也有些类似于分类树了。回归树

当决策树的输出变量（因变量）是分类变量时，叫分类树，而当决策树的输出变量为连续变量时称为回归树。虽然回归树的因变量是连续变量，但叶节点数目是有穷的，因此输出的值也是在这个叶节点上的观测值的平均。回归树不用假定经典回归中的诸如独立性、正态性、线性或者光滑性等等。它对于数量变量和定性变量都同样适用。然而它需要更多的数据来保证合理的结果。例(数据B1.txt).

这是波士顿郊区的房价数据的一部分[1]。原数据有506个街区（观测）及14个变量。我们这里只取了3个变量：人均犯罪率(crim)、较低地位人的比率(lstat)及平均房价(medv，单位：千美元)。我们试图用crim和lstat来预测平均房价medv。首先，我们对变量crim进行对数变换[1]，图10.3是变量crim在变换前后和其它两个变量做散点图的结果。从图中可以看出，变换改变了点的分布极端不均匀的现象。对例9.3数据中的变量crim进行对数变换前后和另外两个变量所做的散点图。

B1=read.table("f:/hepbook/data/Boston.txt",header=T)数据B1.txtlibrary(tree)(b=tree(medv~lstat+crim,B1))Devianceplot(b);text(b)plot(B1$lstat,B1$crim,xlab="lstat",ylab="crim");partition.tree(b,add=TRUE,cex=1.5)思考一下：回归树和分类树有什么异同？回归树和经典回归有什么区别？回归树模型有没有线性、非线性之分，为什么？组合方法：

adaboost、bagging和随机森林

为什么组合？

考虑某人欲竞选某地领导，假定该地有49%的人不支持他。那么，每随机问一个人，都有49%的可能不选他（相当于Bernoulli试验）。如果从该地随机选择1000人来投票，按照简单多数当选的原则，那么他不被选上的概率是多少呢？假定这次投票中不选他的票数服从参数为1000和0.49的二项分布；容易计算，这1000人中有超过半数的人（至少501人）不选他的概率为约0.2532498，远小于某一个人不选举他的概率0.49。n=10000;p=.51s=seq(1,n,2)x=NULLfor(iins)x=c(x,(pbinom(floor(i/2),i,p,lower=F)))plot(s,x,type="l")n=10000;p=.49s=seq(1,n,2)x1=NULLfor(iins)x1=c(x1,(pbinom(floor(i/2),i,p,lower=F)))plot(s,x1,type="l")仅仅是直观描述……par(mfrow=c(1,2));plot(s,x,type="l");plot(s,x1,type="l")仅仅是直观描述……par(mfrow=c(2,3));a=seq(.01,1,.01);n=c(5,11,51,101,1001,9999)for(iin1:6){plot(a,1-pbinom(floor(n[i]/2),n[i],a),type="l",main=substitute(n==q,list(q=n[i])),ylab=expression(paste(群体投票决策概率,"",p[g])),xlab=expression(paste(个体决策概率,"",p)))abline(h=0.5,lty=3);abline(v=0.5,lty=3);segments(0,0,1,1,lty=2)}仅仅是直观描述……仅仅是直观描述……Adaboost

(adaptiveboosting)假定我们的目的是分类。开始可能用的是一种较弱的分类器（即出错率较高的分类器），然后，随着迭代，不断地通过加权再抽样改进分类器，每一次迭代时都针对前一个分类器对某些观测值的误分缺陷加以修正，通常是在（放回）抽取样本时对那些误分的观测值增加权重（相当于对正确分类的减少权重），这样就形成一个新的分类器进入下一轮迭代；而且在每轮迭代时都对这一轮产生的分类器给出错误率。最终结果由各个阶段的分类器的按照错误率加权（权重目的是惩罚错误率大的分类器）投票产生。这就是所谓的“自适应”。Adaboost

Adaboost的缺点是对于奇异点或离群点比较敏感，但其优点是对于过拟合则不那么敏感。Adaboost与boosting(助推法)有很多相似之处，这些方法衍生出来了一系列方法；各种方法都有其优缺点。IRISlibrary(adabag);library(rpart);set.seed(0)samp<-c(sample(1:50,25),sample(51:100,25),sample(101:150,25))a=adaboost.M1(Species~.,data=iris[samp,],mfinal=15,maxdepth=5)a.pred<-predict.boosting(a,newdata=iris[-samp,]);a.pred[-1]a.predt<-predict.boosting(a,newdata=iris[samp,]);a.predt[-1]IRIS：变量重要性barplot(a$importance)例9.5(数据Vehicle.txt)846个观测值，描述了19个变量的数据[1]。我们希望用其中18个关于汽车的度量（数量变量）来把汽车分类。数据中最后一个变量class为已知的类别（四种：bus、opel、saab、van）。例9.5(数据Vehicle.txt)library(mlbench);data(Vehicle);n<-length(Vehicle[,1])samp<-sample(1:n,n/2);a=adaboost.M1(Class~.,data=Vehicle[samp,],mfinal=25,maxdepth=5)a.pred<-predict.boosting(a,newdata=Vehicle[-samp,]);a.pred[-1]b.pred<-predict.boosting(a,newdata=Vehicle[samp,]);b.pred[-1]barplot(a$importance)思考一下：我们用的软件包adabag在输出中除了结果的分类预测之外还有每次迭代时的分类树。以及这些树对每个观测值的“投票。”虽然adaboost的方法描述比较罗唆，但使用起来和其它分类方法没有多大区别。Bagging(bootstrapaggregating)它利用了自助法(bootstrap)放回抽样。它对训练样本做许多次（比如k次）放回抽样，每次抽取和样本量同样的观测值（但由于是放回抽样，大约有三分之一的观测没有抽到）；于是就有了k个不同的样本。然后，对每个样本生成一个决策树。这样，每个树都对一个新的观测值产生一个预测。如果目的是分类，那么由这些树的分类结果的多数（“投票”）产生bagging的分类；如果目的是回归，则由这些树的结果的平均得到因变量的预测值。例9.6(继续例9.5数据Vehicle.txt).n=length(Vehicle[,1])samp=sample(1:n,n/2)a=bagging(Class~.,data=Vehicle[samp,],mfinal=25,maxdepth=5)a.pred=predict.bagging(a,newdata=Vehicle[-samp,]);a.pred[-1]b.pred=predict.bagging(a,newdata=Vehicle[samp,]);b.pred[-1]barplot(a$importance)思考一下：Bagging的各个子分类器是同等权数投票，而adaboost的各个分类器是根据各个分类器的表现加权投票。想想为什么每次自助法抽样有大约三分之一的观测值不出现在自助法样本中。随机森林

随机森林也是进行许多次自助法放回抽样；它的样本数目要大大多于bagging此外，在生成树的时候，在每个节点都仅仅在随机选出的少数变量中选择。不但样本是随机的，每棵树、每个节点的产生都有很大的随机性。随机森林让每个树尽量增长，而且不进行修剪。随机森林

优点它比它产生以前所有的方法都精确。此外，对于大的数据库，它很有效率。它不惧怕很大的维数，即使是数千变量，它也不必删除变量。它还给出分类中各个变量的重要性。随着森林的增长，它产生一个内部无偏的一般误差的估计。它有一个有效的方法来估计缺失值，同时在很大比例数据缺失时仍然保持精确。在一个关于淋巴瘤的基因芯片数据中，变量个数可以达到4682个，而样本量仅有81个；这种数据在经典统计中根本无法处理，Diaconis&Efron(1983)[1]曾经说过，“统计经验表明，基于19个变量和仅仅155个数据点来拟合模型是不明智的。”但是随机森林可以很好地找到重要的基因[2]

例9.7(继续例9.3数据B1.txt).library(randomForest)B1=read.table("f:/hepbook/data/Boston.txt",header=T)#¸ÃÊý¾ÝÔÚMASSattach(B1)w=randomForest(medv~.,data=B1,importance=TRUE,proximity=TRUE)par(mfrow=c(1,2))for(iin1:2)barplot(t(importance(w))[i,],s=0.5)print(w)随机森林：IRISset.seed(117)w=randomForest(Species~.,data=iris,importance=TRUE,proximity=TRUE)par(mfrow=c(1,2));for(iin1:2)barplot(t(importance(w))[i,],s=0.7)随机森林:IRISpar(mfrow=c(1,2))aa=eigen(w$proximity)[[2]]plot(aa[,1:2],pch=(1:3)[as.numeric(iris$Species)],cex=0.7,main="",xlab="",ylab="")MDSplot(w,iris$Species,palette=rep(1,3),pch=as.numeric(iris$Species),xlab="",ylab="")例9.9(数据airquality.txt)1973年纽约市从五月到九月的空气质量数据[1]。变量有臭氧(Ozone，单位ppb)、阳光辐射(Solar.R,单位lang)、风力(Wind,单位mph)、温度(Temp,单位F)、月(Month)、日(Day)。我们希望把臭氧作为因变量，用气象数据作为自变量通过随机森林进行回归。例9.9(数据airquality.txt)ozone.rf=randomForest(Ozone~.,data=airquality,mtry=3,importance=TRUE,na.action=na.omit)par(mfrow=c(1,2));for(iin1:2)barplot(importance(ozone.rf)[,i])(该数据在randomForest中固有)思考一下：随机森林、bagging、adaboost等方法都是基于决策树的组合方法。但这些方法和传统统计思维有很大的不同。你能够说出你的感觉吗？就你的理解，这些方法有什么优缺点？最近邻方法

最近邻方法(nearestneighboralgorithm)可能是所有算法建模中最简单的方法。它基于训练集对测试集进行分类或回归。每个回归或分类问题都有一些自变量，它们组成一个多维空间。首先在空间中假定一个距离。在连续型自变量的情况，通常都是用欧氏空间。在分类问题中，一个测试集点应该属于离它最近的k个训练集点中多数所属于的类型。在k=1的最简单的情况，那么该点的类型应该和与它最近的一个点相同。在回归中，一个测试集点的因变量的预测值应该等于离它最近的k个训练集点相应的因变量值的平均。至于k的选择，一般都用测试集的交叉验证来进行。例9.10(数据radar.txt).关于雷达是否被电离层反射的具有351个观测值的数据[1]。这里自变量有34个，因变量class为分类变量，它取两个值：一个是的确被电离层反射，标为“好”(good)，否则就是“不好”(bad)。我们在这里试图用最近邻方法来建立分类模型用那34个自变量预测这个分类变量。我们还是随机取约一半数据作为训练集，其余为测试集。例9.10(数据radar.txt).library(kknn);data(ionosphere)n=nrow(ionosphere);set.seed(1);test=sample(1:n,n/2)a=kknn(class~.,ionosphere[-test,],ionosphere[test,])table(ionosphere[test,]$class,a$fit)library(mlbench);data(BreastCancer)

w=BreastCancer例3.3(BreastCancer.txt)这是加州大学伯克利分校数据库[1]提供的关于威斯康星大学麦迪逊分校医院的乳腺癌数据(其应用参见Wolberg&Mangasarian,1990和Zhang,J.,1992)。一共有699个观测值和11个变量。目标变量（因变量）是关于病人是良性(benign)还是恶性(malignant)肿瘤的二分变量class，而协变量包括病人的ID、肿块厚度(Cl.thickness)、细胞大小(Cell.size)、细胞形状(Cell.shape)、边缘粘连(Marg.adhesion)、单独上皮细胞大小(Epith.c.size)、裸细胞核(Bare.nuclei)、淡染色质(Bl.cromatin)、正常细胞核(Normal.nucleoli)、分裂激素(Mitoses)。这些协变量都是数量型变量，除了ID之外，都转换成了1到10的整数。整个数据有16个缺失值。我们把缺失值去掉，随机选择三分之一的数据作为测试集，而剩下的三分之二作为训练集。选择k=5，利用triangularkernel作为核函数进行加权KNN分类，对于测试集的预测结果在下面表中：library(mlbench);data(BreastCancer);w=BreastCancerlibrary(kknn)w=w[-1]w=na.omit(w)m<-dim(w)[1]#m=699-16=683set.seed(10)val<-sample(1:m,size=round(m/3),replace=FALSE,prob=rep(1/m,m))w.learn<-w[-val,];dim(w.learn)#455*10w.valid<-w[val,];dim(w.valid)#228*10w.kknn<-kknn(Class~.,k=5,w.learn,w.valid,distance=1,kernel="triangular")summary(w.kknn)fit<-fitted(w.kknn)#228table(w.valid$Class,fit)library(ElemStatLearn)str(mixture.example)#class(mixture.example);attributes(mixture.example)$names[1]"x""y""xnew""prob""marginal""px1"[7]"px2""means"x<-mixture.example$xg<-mixture.example$yx.mod<-lm(g~x)#Figure2.1:plot(x,col=ifelse(g==1,"red","green"),xlab="x1",ylab="x2")coef(x.mod)abline((0.5-coef(x.mod)[1])/coef(x.mod)[3],-coef(x.mod)[2]/coef(x.mod)[3])ghat<-ifelse(fitted(x.mod)>0.5,1,0)#?ifelselength(ghat)sum(ghat==g)1-sum(ghat==g)/length(g)#[1]0.27#TrainingmisclassificationrateTheElementsofStatisticalLearning,DataMining,Inference,andPrediction

byTrevorHastie,RobertTibshiraniandJeromeFriedman.xnew<-mixture.example$xnewdim(xnew)colnames(xnew)library(class)mod15<-knn(x,xnew,g,k=15,prob=TRUE)summary(mod15)#Figure2.2:plot(x,col=ifelse(g==1,"red","green"),xlab="x1",ylab="x2")str(mod15)prob<-attr(mod15,"prob")prob<-ifelse(mod15=="1",prob,1-prob)#probisvotingfractionforwinningclass!#Nowitisvotingfractionforred==1px1<-mixture.example$px1px2<-mixture.example$px2prob15<-matrix(prob,length(px1),length(px2))contour(px1,px2,prob15,levels=0.5,labels="",xlab="x1",ylab="x2",main="15-nearestneighbour")#addingthepointstotheplot:points(x,col=ifelse(g==1,"red","green"))ghat15<-ifelse(knn(x,x,k=15,cl=g)=="1",1,0)sum(ghat15==g)#[1]1691-sum(ghat15==g)/length(g)#[1]0.155#Misclassificationrateforknn(,k=15)#Thenwewanttheplotforknnwithk=1:(Figure2.3)mod1<-knn(x,xnew,k=1,cl=g,prob=TRUE)prob<-attr(mod1,"prob")prob<-ifelse(mod1=="1",prob,1-prob)#probnowisvoting#fractionfor"red"prob1<-matrix(prob,length(px1),length(px2))contour(px1,px2,prob1,level=0.5,labels="",xlab="x1",ylab="x2",main="1-nearestneighbour")#Addingthepointstotheplot:points(x,col=ifelse(g==1,"red","green"))#Reproducingfigure2.4,page17ofthebook:#Thedatadonotcontainatestsample,sowemakeone,#usingthedescriptionoftheoraclepage17ofthebook:Thecenters#isinthemeanscomponentofmixture.example,withgreen(0)first,#sored(1).Foratestsampleofsize10000wesimulate#5000observationsofeachclass.library(MASS)set.seed(123)centers<-c(sample(1:10,5000,replace=TRUE),sample(11:20,5000,replace=TRUE))means<-mixture.example$meansmeans<-means[centers,]mix.test<-mvrnorm(10000,c(0,0),0.2*diag(2))mix.test<-mix.test+meanscltest<-c(rep(0,5000),rep(1,5000))

ks<-c(1,3,5,7,9,11,15,17,23,25,35,45,55,83,101,151)#nearestneighbourstotrynks<-length(ks)misclass.train<-numeric(length=nks)misclass.test<-numeric(length=nks)names(misclass.train)<-names(misclass.test)<-ksfor(iinseq(along=ks)){mod.train<-knn(x,x,k=ks[i],cl=g)mod.test<-knn(x,mix.test,k=ks[i],cl=g)misclass.train[i]<-1-sum(mod.train==factor(g))/200misclass.test[i]<-1-sum(mod.test==factor(cltest))/10000}print(cbind(misclass.train,misclass.test))#Usingpackagemclust02#NotethatthispackageisnolongeronCRAN,#butmustbesearchedinthearchives.

if(require(mclust02)){x<-mixture.example$xg<-mixture.example$yxnew<-mixture.example$xnewpx1<-mixture.example$px1px2<-mixture.example$px2mix.mclust<-mclustDA(x,g,xnew,G=1:6,verbose=TRUE)mix.mclust}#endrequire(mclust02)载入需要的程辑包：mclust02Failedwitherror:‘mclust02'不是有效的已经安装了的程序包’失败

#Figure2.4plot(misclass.train,xlab="NumberofNN",ylab="Testerror",type="n",xaxt="n")axis(1,1:length(ks),as.character(ks))lines(misclass.test,type="b",col='blue',pch=20)lines(misclass.train,type="b",col='red',pch=20)legend("bottomright",lty=1,col=c("red","blue"),legend=c("train","test"))#Figure2.5prob<-mixture.example$probprob.bayes<-matrix(prob,length(px1),length(px2))contour(px1,px2,prob.bayes,levels=0.5,labels="",xlab="x1",ylab="x2",main="Bayesdecisionboundary")points(x,col=ifelse(g==1,"red","green"))思考一下：最近邻方法是最简单最直观的方法。似乎俗语“物以类聚，人以群分”可以描述这个方法。最近邻方法最常用作法的不是让最近的k个点简单地投票，而是加权，即让更接近的点的投票分量越重。你可以思考一下加权的方法。人工神经网络

人工神经网络(ArtificialNeuralNetworks)是对自然的神经网络的模仿；它可以有效地解决很复杂的有大量互相相关变量的回归和分类问题。下面就是一个有两个自变量（输入）和一个因变量（输出）的神经网络的示意图。左边代表自变量的两个节点(node)形成输入层(inputlayer)，中间三个节点形成隐藏层(hiddenlayer)，最右边的一个节点属于输出层(outputlayer)，代表因变量。这些节点按照箭头连接。最常用的隐藏层激活函数为S型logistic函数关键在于权重，需要计算误差，以反馈回上一层节点，调整权重复杂的人工神经网络library(nnet);library(mlbench);data(Vehicle);n=length(Vehicle[,1]);set.seed(1);samp=sample(1:n,n/2);b=class.ind(Vehicle$Class);test.cl=function(true,pred){true<-max.col(tr