Algae(海藻)例子.ppt

1,Algae(海藻)例子,algaetrain-read.table(f:/讲课/SJM/Analysis.txt,header=T,dec=.,na.strings=c(XXXXXXX),2,algaetrain-read.table(f:/讲课/SJM/Analysis.txt,header=T,dec=.,na.strings=c(XXXXXXX)algaetest1-read.table(f:/讲课/SJM/Testdata.txt,header=T,dec=.,na.strings=c(XXXXXXX)algaetest2-read.table(f:/讲课/SJM/Solutions.txt,header=T,dec=.,na.strings=c(XXXXXXX)algaetest-cbind(algaetest1,algaetest2)algae-rbind(algaetrain,algaetest)summary(algae)#对algae数据集的基本信息特征的描述,3,par(mfrow=c(1,2)hist(algae$mxPH,prob=T)hist(algae$mxPH,prob=T,xlab=,main=HistogramofMaximumPHvalue,ylim=0:1)lines(density(algae$mxPH,na.rm=T)#在直方图中添加密度线rug(jitter(algae$mxPH)#画出变量mxPH在X轴附近的真实值,4,5,library(lattice)bwplot(sizea1,data=algae,ylab=RiverSize,xlab=AlgaA1)#,6,缺失值的处理（以对训练集algaetrain的缺失值处理为例）,7,algaetrain!complete.cases(algaetrain),#显示algaetrain数据集中含有缺失值项nrow(algaetrain!complete.cases(algaetrain),)#缺失值项的总行数,8,方法1：直接将含有缺失值项的水样删除,algaetrain1-na.omit(algaetrain),9,方法2：运用常用的均值和中位数的方法进行插补,10,algaetrain2mean-algaetrainfor(rinwhich(!complete.cases(algaetrain)algaetrain2meanr,which(is.na(algaetrainr,)-apply(data.frame(algaetrain,which(is.na(algaetrainr,),2,mean,na.rm=T)#对含有缺失值的记录用均值进行插补,algaetrain2median-algaetrainfor(rinwhich(!complete.cases(algaetrain)algaetrain2medianr,which(is.na(algaetrainr,)-apply(data.frame(algaetrain,which(is.na(algaetrainr,),2,median,na.rm=T)#对含有缺失值的记录用中位数进行插补,11,通过对含缺失值项的变量同其他变量的相关性的大小建立回归方程来插补,12,cor(algaetrain,4:18,use=complete.obs)#变量间的相关性，计算时不考虑缺失值symnum(cor(algaetrain,4:18,use=complete.obs)#简单明显示出变量间的相关,13,14,从变量之间的相关程度图中，我们可以发现绝大部分变量之间相关程度不高，但变量NH4和NO3、变量PO4、oPO4之间相关程度很大，尤其是变量PO4、oPO4之间相关程度达到91.2%。由于数据集algaetrain中第62条和第199条含有的缺失项过多，我们可以考虑将其从数据集中直接删除，而变量NH4和NO3的缺失值集中在数据集algaetrain中第62条和第199条中，故不需要再考虑对其插补。下面来考虑对变量PO4、oPO4建立回归方程来对缺失值进行插补。algaetrain3l-lm(oPO4PO4,data=algaetrain)#对变量PO4、oPO4建立回归方程summary(algaetrain3l)#建立函数fillPO4，对缺失项进行回归插补,15,fillPO4-function(oP)if(is.na(oP)return(NA)elsereturn(oP+15.6142)/0.6466)algaetrain3-algaetrainalgaetrain3is.na(algaetrain$PO4),PO4-sapply(algaetrainis.na(algaetrain$PO4),oPO4,fillPO4)#对algaetrain中的PO4缺失值项进行回归插补,16,方法4：从数据集中选取若干条水样成分和它相似的水样(即和它在空间欧式距离最短的n条水样)，求其中位数进行插补,17,#由于数据集algaetrain中第62条和第199条中含有过多的缺失值项，我们可#以考虑直接将这两条水样删除algaetraind-algaetrain-c(62,199),algaetrain4-algaetraindlibrary(cluster)#调用R语言中的cluster包dist.mtx-as.matrix(daisy(algaetrain4,stand=T)#计算这198个水样的空间距离，#并表示为矩阵的形式，将结果赋值给变量dist.mtxcentral.value-function(x)if(is.numeric(x)median(x,na.rm=T)elseif(is.factor(x)levels(x)which.max(table(x)elsef-as.factor(x)levels(f)which.max(table(f)#建立函数central.value，如果x是数值型的，则取x的中位数赋值给central.value，其他类型的转化为因子型，再输出因子型的最大的那一项的名称。for(rinwhich(!complete.cases(algaetraind)algaetrain4r,which(is.na(algaetraindr,)-apply(data.frame(algaetraindc(eger(names(sort(dist.mtxr,)2:11),which(is.na(algaetraindr,),2,central.value)#将algaetrain4中所含的缺失值项选取在空间上离它最近的2到11个点，取#其距离为central.value那个值来对其列变量中缺失项赋值。#注意，这里用了data.frame，不用等于直接对数据集algae4赋值，否则会#出现变量个数与赋值的个数不一致的报错。,18,模型的建立,#研究藻类a1同哪些因素有关#将以方法（4）处理后的训练数据集algaetrain4作为最终待处理的数据集。并将其赋与新的数据集，命名为：clean.algae。clean.algae-algaetrain4,19,传统统计模型的建立,20,#因为变量season、size、speed为定性变量，在建立回归时，应转化为factor型。lm.a1=lm(a1factor(season)+factor(size)+factor(speed)+mxPH+mnO2+Cl+NO3+NH4+oPO4+PO4+Chla,data=clean.algae,1:12)summary(lm.a1),21,#lm.a1的方差分析anova(lm.a1),22,#从上述方差分析的输出结果中，我们可以看出，因子season对模型拟合优度的贡献最小，可以考虑将这一项从方程中剔除。lm.a2-update(lm.a1,.-factor(season)summary(lm.a2),23,#从lm.a2和lm.a1的AdjustedR-squared比较，0.3281比0.3217略微有所提高，但拟合优度并未得到多大的提高，进一步看看lm.a2和lm.a1的方差比较，anova(lm.a1,lm.a2),#从输出结果，我们可以看出，残差平方和，lm.a2比lm.a1有所下降（-396），但F检验没有通过，说明模型lm.a2和lm.a1差异并不明显。,24,下面用逐步回归法来进行处理（选择方法：全部、前向、后向）,25,lm1both-step(lm.a1,direction=both)lm1forward-step(lm.a1,direction=forward)lm1back-step(lm.a1,direction=backward)final.lm-step(lm.a1)summary(final.lm),26,#下面检验回归方程的残差项是否服从独立同分布#残差项正态性检验及作图（残差图、检验正态性的QQ图等等）shapiro.test(final.lm$residuals),27,par(mfrow=c(1,2)hist(final.lm$residuals)qqnorm(final.lm$residuals),28,从回归方程的输出结果，我们可以看出，AdjustedR-squared仅仅只有0.3333，而且方程中有几项的回归系数并不显著。并且残差项并非正态，说明回归效果不好。进一步我们可以看看它的AIC值。AIC值很大，这说明在这里建立回归方程并不适合。,predictlm-predict(final.lm)cat(逐步回归训练集上的NMSE为:,mean(clean.algae$a1-as.numeric(predictlm)2)/mean(mean(clean.algae$a1)-clean.algae$a1)2),n)#输出逐步回归训练集上的NMSE为:0.6463594,29,#将测试集上的缺失值用上面的方法四来填补，形成最终的测试集algaetestzz=algaetestfor(rinwhich(!complete.cases(algaetest)algaetestzzr,which(is.na(algaetestr,)-apply(data.frame(algaetestc(eger(names(sort(dist.mtxr,)2:11),which(is.na(algaetestr,),2,central.value)nrow(algaetestzz!complete.cases(algaetestzz),),predictlmcs-predict(final.lm,algaetestzz,1:12)cat(逐步回归测试集上的NMSE为:,mean(algaetestzz$a1-as.numeric(predictlmcs)2)/mean(mean(algaetestzz$a1)-algaetestzz$a1)2),n)逐步回归测试集上的NMSE为:0.7041271AIC(final.lm)11703.994,30,用最佳子集（调用leaps包）进行处理,31,library(leaps)lm2leaps-regsubsets(a1factor(season)+factor(size)+factor(speed)+mxPH+mnO2+Cl+NO3+NH4+oPO4+PO4+Chla,data=clean.algae,nbest=1,nvmax=15,really.big=TRUE)#方程中不含常数项共有15个自变量，所以nvmax=15，而nbest为相同个数的自变量的最佳组合排列的方程个数。系统默认为1，即取最好的，因为我要找出所有变量组合中最好的模型，所以为了计算简便，就取系统默认的，因为运算量较大，really.big=TRUE，系统默认really.big=FALSE）。lm2leaps.sum-summary(lm2leaps)names(lm2leaps.sum)1whichrsqrssadjr2cpbicoutmatobj#建立仅有截距项的模型，即规模为0的模型lm2leapsdum-lm(a11,data=clean.algae)lm2leapsdum.sum-summary(lm2leapsdum)names(lm2leapsdum.sum)1calltermsresidualscoefficientsaliasedsigmadf8adj.r.squaredr.squaredcov.unscaledlm2leapsdum.sum$adj10#这说明方程中仅含有常数项的方程比按最佳子集回归得到的15个方程都差，可舍弃#下面我们来看看按最佳子集回归得到的15个方程的拟合效果的比较,32,lm2leaps.sum$rss#15个方程的残差平方和比较par(mfrow=c(2,2)plot(1:15,lm2leaps.sum$rss,type=b,xlab=#subsets,ylab=RSS)lm2leaps.sum$adjr2#15个方程的调整R方比较plot(1:15,lm2leaps.sum$adjr2,type=b,xlab=#subsets,ylab=ADJRSquare)lm2leaps.sum$cp#15个方程的cp比较plot(1:15,lm2leaps.sum$cp,type=b,xlab=#subsets,ylab=CP)lm2leaps.sum$bicplot(1:15,lm2leaps.sum$bic,type=b,xlab=#subsets,ylab=BIC),33,34,#综合上面的信息，我们可以发现子集容量在3以后，调整R方上升的幅度很小；同时，cp值在3以后下降的幅度也很小，而BIC信息准则在3处正好达到最优，也正证实了我的猜想，即当子集容量为3时，方程拟合精度最高。即方程中含四自变量（加常数项）的最好。lm2leaps.sum$which3,35,final.lm2leaps-lm(a1factor(size)+NO3+PO4,data=clean.algae)summary(final.lm2leaps),36,shapiro.test(final.lm2leaps$residuals)AIC(final.lm2leaps),建立回归方程并不合适,37,其他方法,38,从用逐步回归建立的模型的方差膨胀因子或容忍度分析，我们有理由怀疑自变量之间存在着多重共线性的影响，实际上，在前面从变量之间的相关程度图中，我们可以发现绝大部分变量之间相关程度不高，但变量NH4和NO3、变量PO4、oPO4之间相关程度很大，尤其是变量PO4、oPO4之间相关程度达到91.2%。这说明变量间从直观上也可以看出有可能存在着多重共线性的影响。对于多重共线性的处理，最经典的处理方法有四种：principalcomponentsregression,ridgeestimation,shrinkageestimatiorsandpartialleastsquares。下面我将从岭回归、lassoregression两种方法来建模,39,岭回归（MASS）,40,41,library(MASS)ridgelm-lm.ridge(a1factor(season)+factor(size)+factor(speed)+mxPH+mnO2+Cl+NO3+NH4+oPO4+PO4+Chla,data=clean.algae,1:12,lambda=seq(0,200,0.01),model=TRUE)names(ridgelm)1coefscalesInterlambdaymxmGCVkHKBkLWridgelm$lambdawhich.min(ridgelm$GCV)#找到GCV最小时对应的lambda#广义的交叉验证准则GCV，越小越好ridgelm$coef,which.min(ridgelm$GCV)#找到GCV最小时对应的系数,#岭回归方程形式：a1=0.5053926-0.2971772+0.5578304-0.5812274+2.6061114-0.2759249-0.8057665-1.6926636mxPH+1.3573421mnO2-2.2162500Cl-2.4007739NO3+0.3359709NH4-2.3073414oPO4-3.0882571PO4-1.7724788Chla,42,par(mfrow=c(1,2)plot(ridgelm)#画出图形，并作出lambda取59.84时的那条线，以红线表示。abline(v=ridgelm$lambdawhich.min(ridgelm$GCV),col=red)#lamda同GCV之间关系的图形plot(ridgelm$lambda,ridgelm$GCV,type=l)abline(v=ridgelm$lambdawhich.min(ridgelm$GCV),col=green),43,岭回归上训练集上的NMSE方为:0.7320761岭回归上测试集上的NMSE为:0.819326,44,lasso回归（lars）,45,LASSO方法在线性模型中，人们必须选择合适的变量；比如常用的逐步回归法就是选择显著的变量而抛弃那些不显著的。Tibshirani(1996)1提出了一个新的方法来处理变量选择的问题。该方法在模型系数绝对值的和小于某常数的条件下，谋求残差平方和最小。该方法既提供了如子集选择方法那样的可以解释的模型，也具有岭回归那样的稳定性。它不删除变量，但使得一些回归系数收缩、变小，甚至为0。因而，该方法被称为lasso(leastabsoluteshrinkageandselectionoperator，最小绝对值收缩和选择算子2)。,1Tibshirani,R.(1996).“Regressionshrinkageandselectionviathelasso”,J.R.Statist.Soc.B,58(1)267-288.2作为非缩写的英文词，lasso还有套马索的意义，和这里的缩写lasso没有任何关联。把它翻译成“套马索”是不恰当的。,46,47,#因为数据集中含有因子项，所以转化为哑元变量：#对训练集变换t=clean.algaeseason=eger(t,1)size=eger(t,2)speed=eger(t,3)y=as.matrix(cbind(season,size,speed)z=as.matrix(t,4:11)x=cbind(y,z)a=NULL;for(iin1:4)a=cbind(a,(x,1=i)*1)colnames(a)=c(autumn,spring,summer,winter)b=NULL;for(iin1:3)b=cbind(b,(x,2=i)*1)colnames(b)=c(riverlarge,rivermedium,riversmall)c=NULL;for(iin1:3)c=cbind(c,(x,3=i)*1)colnames(c)=c(speedhigh,speedlow,speedmedium)w=as.matrix(clean.algae,4:11)xunlian=cbind(a,b,c,w)#对测试集变换p=as.matrix(algaetestzz,4:11)v=algaetestzzseason=eger(v,1)size=eger(v,2)speed=eger(v,3)m=as.matrix(cbind(season,size,speed)f=as.matrix(v,4:11)n=cbind(m,f)a=NULL;for(iin1:4)a=cbind(a,(n,1=i)*1)colnames(a)=c(autumn,spring,summer,winter)b=NULL;for(iin1:3)b=cbind(b,(n,2=i)*1)colnames(b)=c(riverlarge,rivermedium,riversmall)c=NULL;for(iin1:3)c=cbind(c,(n,3=i)*1)colnames(c)=c(speedhigh,speedlow,speedmedium)ceshi=cbind(a,b,c,p),48,library(lars)attach(clean.algae)lassolm-lars(xunlian,a1)lassolmattributes(lassolm),49,lassolm$Cp,50,plot(lassolm),51,#10折交叉验证并找出有着最小值的CVcv.lar,53,#根据Cp准则得出的lasso得出的方程形式为：a1=5.94855517riversmall-0.03283935Cl-0.11164057NO3-0.04636973PO4-0.02537441Chla#10折交叉验证得出的lasso得出的预测值#10折交叉验证得出的lasso得出的预测值predict11-predict(lassolm,data.matrix(xunlian),mode=fraction,s=best,type=fit)predict1-predict(lassolm,data.matrix(ceshi),mode=fraction,s=best,type=fit)cat(lasso训练集上的NMSE为:,mean(clean.algae$a1-as.numeric(predict11$fit)2)/mean(mean(clean.algae$a1)-clean.algae$a1)2),n)lasso训练集上的NMSE为:0.6473745(依据10折交叉验证)cat(lasso测试集上的NMSE为:,mean(algaetestzz$a1-as.numeric(predict1$fit)2)/mean(mean(algaetestzz$a1)-algaetestzz$a1)2),n)lasso测试集上的NMSE为:0.6969239(依据10折交叉验证),54,#根据Cp准则得出的lasso得出的预测值predict21-predict(lassolm,data.matrix(xunlian),mode=step,s=6,type=fit)cat(lasso训练集上的NMSE为:,mean(clean.algae$a1-as.numeric(predict21$fit)2)/mean(mean(clean.algae$a1)-clean.algae$a1)2),n)lasso训练集上的NMSE为:0.7034299(依据Cp准则)predict2-predict(lassolm,data.matrix(ceshi),mode=step,s=6,type=fit)cat(lasso测试集上的NMSE为:,mean(algaetestzz$a1-as.numeric(predict2$fit)2)/mean(mean(algaetestzz$a1)-algaetestzz$a1)2),n)lasso测试集上的NMSE为:0.7492455依据Cp准则)#所以经过比较，我们可以发现根据10折交叉验证得出的lasso回归方程的预测精度比岭回归和根据Cp准则得出的lasso回归方程都要高。所以在这里可以暂时将其看成是最优的拟合模型。#下面我们来看看上面所建的最好的模型（根据10折交叉验证得出的lasso回归方程）它的残差项的正态性检验。shapiro.test(clean.algae$a1-predict11$fit)#训练集上的残差检验W=0.8999,p-value=2.857e-10shapiro.test(algaetestzz$a1-predict1$fit)#测试集上的残差检验W=0.8612,p-value=3.861e-10,55,从传统统计角度来解决上面问题，也许并不太合适，一方面所建的模型中残差项还存在着许多未被提取的信息，另一方面模型的拟合优度和预测精度都不高。下面我们可以换一种思维方式，绕过传统统计的严格的假定（“黑匣子”），运用数据挖掘的一些方法来进行处理。,56,统计应用领域的一些新方法,57,学习算法的不稳定性是指如果训练集有较小的变化，学习算法产生的预测函数将发生较大变化。Breiman指出，决策树、神经网络是不稳定的，而k最近邻方法是稳定的。所以首先，我们可以先从一个稳定的学习算法K最近邻方法来看看对该数据集的拟合效果和预测精度，再从不稳定的算法如决策树、随机森林以及针对不稳定的算法判定树的两种“提升”Boosting和Bagging算法（Breiman指出，稳定性是Boosting和Bagging能否提高预测准确率的关键因素：它们对不稳定的学习算法能提高预测的准确度，而对稳定的学习算法效果不明显，有时甚至使预测精确度降低。）。,58,K近邻回归（kknn包）,主要用于分类。不用拟合模型，对于某点x，找到k个近邻点，然后根据选票，决定该点之类别。,59,library(kknn)a1x=as.matrix(a1)colnames(a1x)=c(a1)xunlianzz-cbind(xunlian,a1x)a1c=as.matrix(algaetest2$a1)colnames(a1c)=c(a1)ceshizz=cbind(ceshi,a1c)knn1lm-kknn(a1.,as.data.frame(xunlianzz),as.data.frame(ceshizz),k=1,distance=1,kernel=rectangular)attributes(knn1lm)cat(knn1测试集上的预测误差为:,1-mean(knn1lm$fitted.values=algaetest2$a1),n)knn1测试集上的预测误差为:0.8357143knn10lm-kknn(a1.,as.data.frame(xunlianzz),as.data.frame(ceshizz),k=10,distance=1,kernel=rectangular)cat(knn10测试集上的预测误差为:,1-mean(knn10lm$fitted.values=algaetest2$a1),n)knn10测试集上的预测误差为:1,60,可以看出，K近邻算法在此处也许并不适应，因为K近邻算法主要适用于因变量为多维标度变量（比如，季节可分别春、夏、秋、冬四季），当因变量为分组类型的话，K近邻回归的效果往往会很好，但此处，a1为数值变量，如果能将a1分别按某一尺度分为大、中、小三种类型再用K近邻回归，往往效果会更好。,61,回归树,下面我们通过另一类的模型回归树（Regressiontrees）来挖掘藻类a1同哪些因素有关。因为回归树能自己处理缺失值，所以我们可以直接利用原数据集algae（只需要剔出含缺失项过多的第62条和第199条水样）。,62,library(rpart)#调用rpart包algaetree-algaetrain-c(62,199),rt.a1-rpart(a1.,data=algaetree,1:12)rt.a1,63,plot(rt.a1,uniform=T,branch=1,margin=0.1,cex=0.9)#画出回归树text(rt.a1,cex=0.75)#在树中显示分枝的信息。,64,printcp(rt.a1)#显示回归树rt.a1每一步得出的sub-trees的详细信息,65,rt2.a1-prune(rt.a1,cp=0.07)#以cp=0.07为尺度对树rt.a1进行修剪rt2.a1plot(rt2.a1,uniform=T,branch=1,margin=0.1,cex=0.9)#画出修剪后的回归树text(rt2.a1,cex=0.75)#在树中显示分枝的信息。,66,treeprexl=predict(rt2.a1,algaetree)cat(tree训练集上的NMSE为:,mean(algaetree$a1-as.numeric(treeprexl)2)/mean(mean(algaetree$a1)-algaetree$a1)2),n)tree训练集上的NMSE为:0.5223749treeprecs=predict(rt2.a1,algaetest)cat(tree测试集上的NMSE为:,mean(algaetest$a1-as.numeric(treeprecs)2)/mean(mean(algaetest$a1)-algaetest$a1)2),n)tree测试集上的NMSE为:0.5931156,67,Boosting(（调用mboost包）,68,69,library(mboost)#调用mboost包mboostlm=blackboost(a1.,control=boost_control(mstop=50),data=algaetree,1:12)attributes(mboostlm)t=predict(mboostlm)cat(mboost训练集上的NMSE为:,mean(algaetree$a1-as.numeric(t)2)/mean(mean(algaetree$a1)-algaetree$a1)2),n)mboost训练集上的NMSE为:0.3943555t1=predict(mboostlm,algaetest,1:12)cat(mboost测试集上的NMSE为:,mean(algaetest$a1-as.numeric(t1)2)/mean(mean(algaetest$a1)-algaetest$a1)2),n)mboost测试集上的NMSE为:0.5679551,70,Bagging(bootstrapaggregation),71,Bagging是Breiman提出的与Boosting相似的技术。Bagging技术的主要思想是给定一弱学习算法和一训练集11(,),.,(,)nnxyxy。让该学习算法训练多轮，每轮的训练集由从初始的训练集中随机取出的n个训练例组成，初始训练例在某轮训练集中可以出现多次或根本不出现。训练之后可得到一个预测函数序列：1,.,thh，最终的预测函数H对分类问题采用投票方式，对回归问题采用简单平均方法对新示例进行判别。Bagging与Boosting的区别在于Bagging的训练集的选择是随机的，各轮训练集之间相互独立，而Boosting的训练集的选择不是独立的，各轮训练集的选择与前面各轮的学习结果有关；Bagging的各个预测函数没有权重，而Boosting是有权重的；Bagging的各个预测函数可以并行生成，而Boosting的各个预测函数只能顺序生成。对于象神经网络这样极为耗时的学习方法，Bagging可通过并行训练节省大量时间开销。以下我将通过R语言中的ipred包运用Bagging算法对该数据集进行分析研究。,72,library(ipred)bagging.a1=bagging(algaetree$a1.,data=algaetree,1:12,coob=T,control=rpart.control(xval=10)attributes(bagging.a1)baggingtrain=predict(bagging.a1,algaetree,1:12)baggingpre=predict(bagging.a1,algaetest,1:12)baggingprecat(Bagging训练集上的NMSE为:,mean(algaetree$a1-as.numeric(baggingtrain)2)/mean(mean(algaetree$a1)-algaetree$a1)2),n)Bagging训练集上的NMSE为:0.3562546cat(Bagging测试集上的NMSE为:,mean(algaetest$a1-as.numeric(baggingpre)2)/mean(mean(algaetest$a1)-algaetest$a1)2),n)Bagging测试集上的NMSE为:0.4805426bagging.a1$mtrees#大量的树（按照缺省值为25个）,73,随机森林（RandomForest）,74,library(randomForest)#调用randomForest包#因为randomForest不能自动处理缺失值，得插补algaetestt=algaetest#对训练集中含有缺失项的用方法四进行插补for(rinwhich(!complete.cases(algaetest)algaetesttr,which(is.na(algaetestr,)-apply(data.frame(algaetestc(eger(na