课件05fz多元数据分析

多元数据分析经典=>可视化实数轴上的数据:经典多元分析内容回顾本书的所有内容都是关于多元数据的,这里的标题主要是将就经典``多元分析''的名称和内容.经典的多元分析是以对数据的正态性假定为基础的.包括聚类分析,主成份分析及因子分析,典型相关分析,判别分析.有的教科书也包括不那么经典的对应分析.经典:主成份分析及因子分析主成份分析及因子分析的目的都是降维,就是把相关的变量数目减少,用较少的变量来取代原先变量,而这些新变量为原始变量的线性组合.如果原始变量互相正交,那么这两种方法都不可用.这两个方法的区别不在于数据或研究目的,而是数学模型的区别.在数学意义上,主成份分析是因子分析的特例.一般来说,因子分析的效果有可能优于主成份分析.主成份分析优点是原理简单,计算容易.因子分析实际上是后面要引进的图模型或隐变量模型的特例.由于很强的关于正态性及协方差矩阵的假定,因子分析甚至可以只要样本协方差矩阵而非原始数据也可以运行,这与用协方差解图模型的方法一致.当然,较强的假定往往使得模型距离现实更远,也更有理由对结果作出质疑.在做主成份分析及因子分析时,最主观的,也是对结果影响最大的就是参与建模的变量选择,这方面怎么强调也不过分.例5.1洛杉矶街区(LA.Neighborhoods.csv).}该数据为美国普查局2000年数据.一共有110个街区,15个变量,变量情况如下表:主成份分析w=read.csv("LA.Neighborhoods.csv")#读入数据w$density=w$Population/w$Area#增加人口密度变量u=w[,-c(12:15)]#去掉人口,面积,经纬度变量a=eigen(cor(scale(u[-1])))#解标准化后的数据的特征值-特征向量问题(cca=(a$va)/sum(a$va))#贡献率(ca=cumsum(a$va)/sum(a$va))#累计贡献率par(mfrow=c(1,2))#画贡献率图(崖底碎石图-ScreePlot)及累积贡献率图:plot(1:11,a$va,type="o",pch=17,col=4,main="ScreePlot",xlab="ComponentNumber",ylab="EigenValue")plot(1:11,ca,type="o",pch=17,col=4,main="CumulativeContribution",xlab="ComponentNumber",ylab="CumulativeContribution")par(mfrow=c(1,1))(loadings=sweep(a$vec,2,sqrt(a$value),"*"))->b#计算载荷par(mfrow=c(1,2))plot(b[,1:2],type="n",xlab="Component1(42%)",ylab="Component2(17%)",xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),main="Loadings")text(b[,1],b[,2],names(u[,-1]))abline(h=0);abline(v=0)plot(b[,3:4],type="n",xlab="Component3(11%)",ylab="Component4(9%)",xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),main="Loadings")text(b[,3],b[,4],names(u[,-1]))abline(h=0);abline(v=0)par(mfrow=c(1,1))par(mfrow=c(1,2))sc=as.matrix(scale(u[,-1]))%*%a$ve[,1:4]#计算得分plot(sc[,1],sc[,2],type="n",ylim=c(-6,6),xlim=c(-7,6),main="SamplePrincipalComponents",xlab="Component1",ylab="Component2")text(sc[,1],sc[,2],u[,1],cex=.4);abline(v=0,col=2);abline(h=0,col=2)plot(sc[,3],sc[,4],type="n",ylim=c(-6,6),xlim=c(-7,6),main="SamplePrincipalComponents",xlab="Component3",ylab="Component4")text(sc[,3],sc[,4],u[,1],cex=.4);abline(v=0,col=2);abline(h=0,col=2)par(mfrow=c(1,1))因子分析#下面因子分析时默认旋转rotation="varimax":(a=factanal(factors=2,scale(u[,-1]),scores="regression"))(a$loadings);(a$scores)par(mfrow=c(1,2))plot(a$loadings[,1:2],type="n",ylim=c(-1.1,1.1),xlim=c(-1.1,1.2),xlab="Factor1",ylab="Factor2",main="Loadings")abline(h=0);abline(v=0)text(a$loadings[,1],a$loadings[,2],labels=s(a$loadings),cex=.6)plot(a$loadings[,3:4],type="n",ylim=c(-1.1,1.1),xlim=c(-1.1,1.2),xlab="Factor3",ylab="Factor4",main="Loadings")abline(h=0);abline(v=0)text(a$loadings[,3],a$loadings[,4],labels=s(a$loadings),cex=.6)plot(a$scores[,1:2],type="n",ylim=c(-2,2.2),xlim=c(-3.5,3),xlab="Factor1",ylab="Factor2",main="FactorScores")abline(h=0);abline(v=0)text(a$scores[,1],a$scores[,2],labels=u[,1],cex=.4)plot(a$scores[,3:4],type="n",ylim=c(-2,4),xlim=c(-5,1.5),xlab="Factor3",ylab="Factor4",main="FactorScores",)abline(h=0);abline(v=0)text(a$scores[,3],a$scores[,4],labels=u[,1],cex=.4)分层聚类及k均值聚类分层聚类w=read.csv("LA.Neighborhoods.csv")w=data.frame(w,density=w$Population/w$Area)u=w[,c(1,2,5,6,11,16)]#选择变量hh=hclust(dist(scale(u[,-1])),"complete")#对标准化的数据做分层聚类,聚类方法选的"complete"plot(hh,labels=u[,1],cex=.6)#画树状图id=identify(hh)#手工分成5份K均值聚类a=kmeans(scale(u[,-1]),5);ppp=c(7,17,19,21)plot(w[a$cluster==1,14:15],pch=1,col=1,xlim=c(-118.7,-118.2),ylim=c(33.73,34.32),main="LosAngeles")for(iin2:5)points(w[a$cluster==i,14:15],pch=ppp[i-1])legend("bottomleft",pch=c(1,ppp),paste("Cluster",1:4))典型相关分析data(oliveoil,package="pls");attach(oliveoil)#数据library(CCA)(res.cc=cc(sensory,chemical))$xcoef[,1][,2][,3][,4][,5]yellow-0.2471010.38900.303910.101170.06152green-0.1772120.34450.213460.061640.03562brown0.1011210.2915-0.064650.231420.06991glossy0.0369670.1474-0.469080.165790.60202transp-0.005204-0.19000.27997-0.10092-0.44396syrup-0.066060-0.43470.43096-0.179620.11612$ycoef[,1][,2][,3][,4][,5]Acidity2.28933-6.83138-1.8312-8.01461.5270Peroxide0.06877-0.087130.26350.14130.5405K2324.31628-6.11597-1.6677-5.1552-8.8541K270-2.4705183.26295-16.069841.004725.1343DK-400.39944636.94814344.9703152.6203-64.4116$scores$corr.X.xscores#相当于cor(sensory,res.cc$scores$xscores)[,1][,2][,3][,4][,5]yellow-0.5523-0.349110.028000.7169-0.12568green0.44160.40651-0.03411-0.74390.13815brown0.81020.095320.359920.44150.09558glossy-0.5475-0.22618-0.480030.40140.04302transp-0.4898-0.24276-0.400150.4880-0.18182syrup0.6588-0.025620.50117-0.25230.47694$scores$corr.Y.xscores#相当于cor(chemical,res.cc$scores$xscores)[,1][,2][,3][,4][,5]Acidity0.027490.26653-0.08334-0.5063530.09375Peroxide0.702260.090850.55018-0.0056380.04450K2320.778040.218400.38155-0.084279-0.07060K2700.729140.50102-0.08460-0.1430350.03332DK-0.012540.372790.56679-0.321762-0.03308$scores$corr.X.yscores#相当于cor(sensory,res.cc$scores$yscores)[,1][,2][,3][,4][,5]yellow-0.5393-0.293150.023050.4109-0.03593green0.43120.34136-0.02807-0.42630.03949brown0.79120.080040.296260.25300.02732glossy-0.5346-0.18993-0.395120.23010.01230transp-0.4782-0.20385-0.329370.2797-0.05197syrup0.6433-0.021510.41253-0.14460.13634$scores$corr.Y.yscores#cor(chemical,res.cc$scores$yscores)[,1][,2][,3][,4][,5]Acidity0.028150.3174-0.1013-0.8835380.3280Peroxide0.719170.10820.6684-0.0098370.1557K2320.796780.26010.4635-0.147060-0.2470K2700.746710.5966-0.1028-0.2495830.1166DK-0.012840.44390.6886-0.561445-0.1157对应分析library(MASS);v=caith#数据来源colnames(v)=paste(colnames(v),"hair")rownames(v)=paste(rownames(v),"eye")(cc=corresp(v,nf=2))#对应分析及结果输出biplot(cc,xlim=c(-.55,.75))#画图非经典多元数据分析:可视化主成份分析w=read.csv("LA.Neighborhoods.csv")u=w[,-c(1,12:15)]library(FactoMineR)res.pca<-PCA(u,quanti.sup=7:10)#补充变量(测试)为7:10个对应分析library(FactoMineR)data(glass,package="anacor")names(glass)=tolower(names(glass))res.ca<-CA(glass)res.ca<-CA(glass,col.sup=5,row.sup=6)多重对应分析例5.4汽车试验(car.test.csv).}这是软件包\verb+rpart+的汽车试验数据,有60个观测,8个变量:价格(Price):单位美元,国家(Country):产地名称,可靠性(Reliability):1-5评分,每加仑英里数(Mileage),类型(Type):6种类型,重量(Weight):单位磅,发动机容量(Disp):升,马力(HP).该数据由Johnetal(1992)引自{\itConsumerReports},April,1990,pp.235–288.data(car.test.frame,package="rpart")library(FactoMineR)#选择两个补充变量,一个数量变量,另一个分类变量:res.mca=MCA(car.test.frame,quanti.sup=7,quali.sup=2)plot(res.mca,invisible=c("var","quali.sup","quanti.sup"),cex=0.7)#选择不显示变量plotellipses(res.mca,keepvar=c(2,5))#画置信椭圆多重因子分析多重因子分析(multipleFactorAnalysis)意在处理分群结构的数据.这些群都由不同类型的变量组成,可以是数量变量,也可以是分类变量.它们都需要同时进行分析.这种数据会出现在许多领域,比如基因领域,传感和物理化学领域,代码的比较等等.对于结构的考虑使得每一群变量的影响得到平衡,考虑了各组变量之间的联系,既给出了传统的图形,也标出了其特殊性.例5.5葡萄酒(wine.csv).有21个观测,31个变量.第一个变量为酒的(三种)来源标签,第二列为土壤(四个水平),其他都是理化指标.library(FactoMineR)data(wine)aa=MFA(wine,group=c(2,5,3,10,9,2),type=c("n",rep("s",5)),ncp=5,name.group=c("orig","olf","vis","olfag","gust","ens"),num.group.sup=c(1,6),graph=FALSE)#这里第1和第6群为补充群plot(aa,choix="ind")#图5.12左plot(aa,choix="ind",partial="all")#图5.12右plot(aa,choix="group")#图5.13左plot(aa,choix="var",habillage="group")#图5.13右plot(aa,choix="axes")#图5.14左plot(1:nrow(aa$eig),aa$eig[,1],type="o",pch=17,col=4,main="Eigenvalues",xlab="DimensionNumber",ylab="Eigenvalues")#图5.14右分层多重因子分析hierar<-list(c(2,5,3,10,9,2),c(4,2))res.hmfa<-HMFA(wine,H=hierar,type=c("n",rep("s",5)))基于主成份分析的聚类例5.6生产(produc.csv).}这是美国48个州从1970到1986年17年的经济数据,变量包括州名(state),年(year),私人资本(pcap),道路(hwy),上下水系统(water),其他公共设施(util),公共资本(pc),州生产总值(gsp),非农村就业(emp),失业率(unemp).library(FactoMineR)data(Produc,package="Ecdat")aa=unique(Produc[,1])z=NULL;for(iin1:48)z=rbind(z,apply(Produc[Produc[,1]==aa[i],-(1:2)],2,mean))s(z)=aares.pca=PCA(z,ncp=5,scale.unit=TRUE,graph=FALSE)#主成份分析res.hcpc=HCPC(res.pca,nb.clust=-1,conso=0,min=3,max=10)#上面nb.clust=-1意味着自动选聚类数目(0是手工选)多元数据的关联规则分析例5.7收入(eESL.csv).}这是软件包\verb+arules+\footnote{MichaelHahsler,BettinaGruenandKurtHornik(2011).arules:MiningAssociationRulesandFrequentItemsets.Rpackageversion1.0-7.}所带数据,来自Hastieetal(2001),有8993个观测值及14个变量.该数据抽自普查数据,变量为收入(e):分为从$[0,10)$到$[75,\infty)$(千美金)的每5千美元一个区间的9个区间,性别(sex):男(male),女(female),婚姻状态(maritalstatus):已婚(married),同居(cohabitation),离婚(divorced),寡居(widowed),单身(single),年龄段(age):14-17,18-24,...,以后10年一段到65+,教育(education):不到9年级9-11年纪,中学毕业(highschoolgraduate),大学肄业(college(1-3years)),大学毕业(collegegraduate),研究生(graduatestudy),职业(occupation):专业人员/经理(professional/managerial),销售(sales),劳力(laborer),职员(clerical/service),家务(homemaker),学生(student),军人(military),退休(retired),失业(unemployed),在旧金山湾区年限(yearsinbayarea):$1,\1-3\4-6\7-10\>10,$未婚者双份收入(duales):是(yes),否(no),家中人口(numberinhousehold):从1到9以上,孩子个数(numberofchildren):从1到9个以上,家庭状况(householderstatus):自有房子(own),租房(rent),和父母或家庭住(livewithparents/family),住房类型(typeofhome):单独房子(house),公寓(condominium),单元房(apartment),移动房屋(mobileHome),其他(other),宗族(ethnicclassification):印第安人(americanindian),亚裔(asian),黑人(black),印度裔(eastindian),美洲西裔(hispanic),太平洋岛民(pacificislander),白人(white),其他(other),母语(languageinhome):英语(english),西班牙语(spanish),其他(other).这里用软件包\verb+arules+来做关联规则分析,首先对数据做整理,下面的代码先把缺失数据删除剩下6876个.由于各个分类变量水平太多,把它们做如下简化:把收入从9个水平换为两个水平(0-40000及40000以上),把年龄换成两个水平(14-35及35以上),把教育换成两个水平(没大学毕业及大学毕业),把在旧金山湾区的年份换成1-9年和9年以上,把家庭人口换成1个及多于1个,把孩子数目换成没有孩子及有孩子.注意,关联规则处理的都是二分变量,如果是连续变量,则应该换成离散型的,而如果不是二分变量,则软件会把它们变成二分变量,比如把有三个水平的变量,自动换成三个以水平为主的两分变量(每个都有``是''与``不是''两个水平).library(arules)data(eESL)data("eESL")eESL<-eESL[complete.cases(eESL),]#不要缺失值

eESL[["e"]]=factor((as.numeric(eESL[["e"]])>6)+1,levels=1:2,labels=c("$0-$40,000","$40,000+"))eESL[["age"]]<-factor((as.numeric(eESL[["age"]])>3)+1,levels=1:2,labels=c("14-34","35+"))eESL[["education"]]=factor((as.numeric(eESL[["education"]])>4)+1,levels=1:2,labels=c("nocollegegraduate","collegegraduate"))eESL[["yearsinbayarea"]]<-factor((as.numeric(eESL[["yearsinbayarea"]])>4)+1,levels=1:2,labels=c("1-9","10+"))eESL[["numberinhousehold"]]<-factor((as.numeric(eESL[["numberinhousehold"]])>3)+1,levels=1:2,labels=c("1","2+"))eESL[["numberofchildren"]]<-factor((as.numeric(eESL[["numberofchildren"]])>1)+0,levels=0:1,labels=c("0","1+"))e<-as(eESL,"transactions")summary(e)itemFrequencyPlot(e,support=0.05,s=0.8)#图5.18fsets=eclat(e,parameter=list(support=0.05,maxlen=10))#频繁项集inspect(sort(fsets,by="support")[1:10])itemssupport1{languageinhome=english}0.91288542{education=nocollegegraduate}0.70520653{numberinhousehold=1}0.69182664{ethnicclassification=white}0.66972085{ethnicclassification=white,languageinhome=english}0.65954046{numberinhousehold=1,languageinhome=english}0.64950557{yearsinbayarea=10+}0.64659698{education=nocollegegraduate,languageinhome=english}0.63438059{e=$0-$40,000}0.622454910{numberofchildren=0}0.6218732rules=apriori(e,parameter=list(support=0.01,confidence=0.01))x=subset(rules,subset=rhs%in%"e=$40,000+"&lift>1.2)inspect(sort(x,by="lift")[1:3])关联规则的可视化library(arulesViz)data(Groceries)rules<-apriori(Groceries,parameter=list(support=0.001,confidence=0.5))rules##scatterplotplot(rules)##try:sel<-plot(rules,interactive=T