05FZ多元數(shù)據(jù)分析

多元數(shù)據(jù)分析=>可視化1實數(shù)軸上的數(shù)據(jù):經(jīng)典多元分析內(nèi)容回顧本書的所有內(nèi)容都是關于多元數(shù)據(jù)的,這里的標題主要是將就經(jīng)典``多元分析''的名稱和內(nèi)容.經(jīng)典的多元分析是以對數(shù)據(jù)的正態(tài)性假定為基礎的.包括聚類分析,主成份分析及因子分析,典型相關分析,判別分析.有的教科書也包括不那么經(jīng)典的對應分析.2經(jīng)典:主成份分析及因子分析主成份分析及因子分析的目的都是降維,就是把相關的變量數(shù)目減少,用較少的變量來取代原先變量,而這些新變量為原始變量的線性組合.如果原始變量互相正交,那么這兩種方法都不可用.這兩個方法的區(qū)別不在于數(shù)據(jù)或研究目的,而是數(shù)學模型的區(qū)別.在數(shù)學意義上,主成份分析是因子分析的特例.一般來說,因子分析的效果有可能優(yōu)于主成份分析.主成份分析優(yōu)點是原理簡單,計算容易.因子分析實際上是后面要引進的圖模型或隱變量模型的特例.由于很強的關于正態(tài)性及協(xié)方差矩陣的假定,因子分析甚至可以只要樣本協(xié)方差矩陣而非原始數(shù)據(jù)也可以運行,這與用協(xié)方差解圖模型的方法一致.當然,較強的假定往往使得模型距離現(xiàn)實更遠,也更有理由對結果作出質(zhì)疑.在做主成份分析及因子分析時,最主觀的,也是對結果影響最大的就是參與建模的變量選擇,這方面怎么強調(diào)也不過分.3例5.1洛杉磯街區(qū)(LA.Neighborhoods.csv).}該數(shù)據(jù)為美國普查局2000年數(shù)據(jù).一共有110個街區(qū),15個變量,變量情況如下表:4主成份分析w=read.csv("LA.Neighborhoods.csv")#讀入數(shù)據(jù)w$density=w$Population/w$Area#增加人口密度變量u=w[,-c(12:15)]#去掉人口,面積,經(jīng)緯度變量a=eigen(cor(scale(u[-1])))#解標準化后的數(shù)據(jù)的特征值-特征向量問題(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))5(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))67par(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))89因子分析10#下面因子分析時默認旋轉(zhuǎn)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)11分層聚類及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")#對標準化的數(shù)據(jù)做分層聚類,聚類方法選的"complete"plot(hh,labels=u[,1],cex=.6)#畫樹狀圖id=identify(hh)#手工分成5份1213K均值聚類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))14典型相關分析15data(oliveoil,package="pls");attach(oliveoil)#數(shù)據(jù)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.411616$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.115717對應分析library(MASS);v=caith#數(shù)據(jù)來源colnames(v)=paste(colnames(v),"hair")rownames(v)=paste(rownames(v),"eye")(cc=corresp(v,nf=2))#對應分析及結果輸出biplot(cc,xlim=c(-.55,.75))#畫圖1819非經(jīng)典多元數(shù)據(jù)分析:可視化20主成份分析w=read.csv("LA.Neighborhoods.csv")u=w[,-c(1,12:15)]library(FactoMineR)res.pca<-PCA(u,quanti.sup=7:10)#補充變量(測試)為7:10個21對應分析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)2223多重對應分析例5.4汽車試驗(car.test.csv).}這是軟件包\verb+rpart+的汽車試驗數(shù)據(jù),有60個觀測,8個變量:價格(Price):單位美元,國家(Country):產(chǎn)地名稱,可靠性(Reliability):1-5評分,每加侖英里數(shù)(Mileage),類型(Type):6種類型,重量(Weight):單位磅,發(fā)動機容量(Disp):升,馬力(HP).該數(shù)據(jù)由Johnetal(1992)引自{\itConsumerReports},April,1990,pp.235–288.data(car.test.frame,package="rpart")library(FactoMineR)#選擇兩個補充變量,一個數(shù)量變量,另一個分類變量: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))#畫置信橢圓2425多重因子分析多重因子分析(multipleFactorAnalysis)意在處理分群結構的數(shù)據(jù).這些群都由不同類型的變量組成,可以是數(shù)量變量,也可以是分類變量.它們都需要同時進行分析.這種數(shù)據(jù)會出現(xiàn)在許多領域,比如基因領域,傳感和物理化學領域,代碼的比較等等.對于結構的考慮使得每一群變量的影響得到平衡,考慮了各組變量之間的聯(lián)系,既給出了傳統(tǒng)的圖形,也標出了其特殊性.例5.5葡萄酒(wine.csv).有21個觀測,31個變量.第一個變量為酒的(三種)來源標簽,第二列為土壤(四個水平),其他都是理化指標.26library(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右272829分層多重因子分析hierar<-list(c(2,5,3,10,9,2),c(4,2))res.hmfa<-HMFA(wine,H=hierar,type=c("n",rep("s",5)))30基于主成份分析的聚類例5.6生產(chǎn)(produc.csv).}這是美國48個州從1970到1986年17年的經(jīng)濟數(shù)據(jù),變量包括州名(state),年(year),私人資本(pcap),道路(hwy),上下水系統(tǒng)(water),其他公共設施(util),公共資本(pc),州生產(chǎn)總值(gsp),非農(nóng)村就業(yè)(emp),失業(yè)率(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意味著自動選聚類數(shù)目(0是手工選)313233多元數(shù)據(jù)的關聯(lián)規(guī)則分析例5.7收入(IncomeESL.csv).}這是軟件包\verb+arules+\footnote{MichaelHahsler,BettinaGruenandKurtHornik(2011).arules:MiningAssociationRulesandFrequentItemsets.Rpackageversion1.0-7.}所帶數(shù)據(jù),來自Hastieetal(2001),有8993個觀測值及14個變量.該數(shù)據(jù)抽自普查數(shù)據(jù),變量為收入(income):分為從$[0,10)$到$[75,\infty)$(千美金)的每5千美元一個區(qū)間的9個區(qū)間,性別(sex):男(male),女(female),婚姻狀態(tài)(maritalstatus):已婚(married),同居(cohabitation),離婚(divorced),寡居(widowed),單身(single),年齡段(age):14-17,18-24,...,以后10年一段到65+,教育(education):不到9年級9-11年紀,中學畢業(yè)(highschoolgraduate),大學肄業(yè)(college(1-3years)),大學畢業(yè)(collegegraduate),研究生(graduatestudy),職業(yè)(occupation):專業(yè)人員/經(jīng)理(professional/managerial),銷售(sales),勞力(laborer),職員(clerical/service),家務(homemaker),學生(student),軍人(military),退休(retired),失業(yè)(unemployed),在舊金山灣區(qū)年限(yearsinbayarea):$1,\1-3\4-6\7-10\>10,$未婚者雙份收入(dualincomes):是(yes),否(no),家中人口(numberinhousehold):從1到9以上,孩子個數(shù)(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).34這里用軟件包\verb+arules+來做關聯(lián)規(guī)則分析,首先對數(shù)據(jù)做整理,下面的代碼先把缺失數(shù)據(jù)刪除剩下6876個.由于各個分類變量水平太多,把它們做如下簡化:把收入從9個水平換為兩個水平(0-40000及40000以上),把年齡換成兩個水平(14-35及35以上),把教育換成兩個水平(沒大學畢業(yè)及大學畢業(yè)),把在舊金山灣區(qū)的年份換成1-9年和9年以上,把家庭人口換成1個及多于1個,把孩子數(shù)目換成沒有孩子及有孩子.注意,關聯(lián)規(guī)則處理的都是二分變量,如果是連續(xù)變量,則應該換成離散型的,而如果不是二分變量,則軟件會把它們變成二分變量,比如把有三個水平的變量,自動換成三個以水平為主的兩分變量(每個都有``是''與``不是''兩個水平).library(arules)data(IncomeESL)data("IncomeESL")IncomeESL<-IncomeESL[complete.cases(IncomeESL),]#不要缺失值IncomeESL[["income"]]=factor((as.numeric(IncomeESL[["income"]])>6)+1,levels=1:2,labels=c("$0-$40,000","$40,000+"))IncomeESL[["age"]]<-factor((as.numeric(IncomeESL[["age"]])>3)+1,levels=1:2,labels=c("14-34","35+"))IncomeESL[["education"]]=factor((as.numeric(IncomeESL[["education"]])>4)+1,levels=1:2,labels=c("nocollegegraduate","collegegraduate"))IncomeESL[["yearsinbayarea"]]<-factor((as.numeric(IncomeESL[["yearsinbayarea"]])>4)+1,levels=1:2,labels=c("1-9","10+"))IncomeESL[["numberinh