函數(shù)庫
001 002 003 004 005 006 007 008 009 010 011 012 013 014 015 016 017 018 019 020 021 022 023 024 025 026 027 028 029 030 031 032 033 034 035 036 037 038 039 040 041 042 043 044 045 046 047 048 049 050 051 052 053 054 055 056 057 058 059 060 061 062 063 064 065 066 067 068 069 070 071 072 073 074 075 076 077 078 079 080 081 082 083 084 085 086 087 088 089 090 091 092 093 094 095 096 097 098 099 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 | #*************因子分析-R語言實(shí)現(xiàn)，函數(shù)庫文件**************##****作者：oldlee11***************************************##****版本：v0.1*******************************************##****時(shí)間：2013-1-17*************************************##功能目標(biāo)：原始數(shù)據(jù)變量x1,x2,x3....xn(全體記為X)。通過樣本可以知道某些變量之間有相關(guān)性。# 則計(jì)算出新變量/因子f1,f2,f3....fm（m<n）（全體記為F），這新變量/因子F可以最大程度的表達(dá)原變量X# 由于新變量的個(gè)數(shù)m小于原始變量個(gè)數(shù)n，即降維了。#原理：AF+e=X# A叫因子載荷（loading）。意義：fi(某1個(gè)因子)和xi(某一個(gè)原變量)的相關(guān)系數(shù)，接近1表示fi和xi相關(guān)性強(qiáng):aij=cov(xi,fj)# e叫特殊因子#其它術(shù)語變量：# 公因子方差：F（所有因子）解釋xi（某一個(gè)原始變量）的方差百分比（貢獻(xiàn)）# 特征值：fi（某一個(gè)因子）解釋X（所有原始變量）的方差百分比（貢獻(xiàn)）# 因子得分：在計(jì)算得出了A后，計(jì)算F內(nèi)的樣本數(shù)據(jù)。# 旋轉(zhuǎn)：對(duì)因子載荷進(jìn)行旋轉(zhuǎn)，之后因子載荷各項(xiàng)大的越來越大，小的越來越小，便于劃分因子和原始變量的關(guān)系。#****函數(shù)：factor()#****概要：因子分析法#****輸入：# 名稱 | 數(shù)據(jù)格式# data_frame | 欲分析的數(shù)據(jù) ，數(shù)據(jù)框格式，最好帶名稱，不可以有factor分類數(shù)據(jù)# factors | 欲產(chǎn)生的因子個(gè)數(shù)小于原始數(shù)據(jù)的變量個(gè)數(shù)# scores | 是否進(jìn)行因子得分# rotation | 是否進(jìn)行旋轉(zhuǎn)# loadings.abs.std | 用于從因子載荷中挑出那些原始變量和那些因子相關(guān)的標(biāo)準(zhǔn)（相對(duì)系數(shù)的絕對(duì)值的最低標(biāo)準(zhǔn)）#****輸出：# factanal(..)產(chǎn)生的結(jié)果# factanal(..)$score并非數(shù)據(jù)框格式，需要as.data.frame(factanal(..)$score)轉(zhuǎn)化一下才是數(shù)據(jù)框格式。factor<-function(data_frame,factors=2,scores="none",rotation="none",loadings.abs.std=0.6){ sol<-factanal(~.,data=data_frame,factors=factors,scores=scores,rotation=rotation)#使用最大似然法進(jìn)行的因子分析。 print("==========================") print("==== 因子分析結(jié)果如下 ====") print("==========================") print("模型：X=AF+e") print("======================================================================") print("1 A即因子載荷loadings:") print(" 每個(gè)數(shù)據(jù)(aij)表示了原始變量xi和因子變量fj的相關(guān)系數(shù)cov.值約接近+1或-1，約相關(guān)") for(i in 1:factors){ print(paste(" 因子",i,"同原始變量的相關(guān)性系數(shù)")) print(rev(sort(loadings(sol)[,i]))) } #### 給出因子和原始變量的可能關(guān)系#### print(" 您可以通過以上數(shù)據(jù)查看接近+1或者-1的數(shù)據(jù)，以說明某一因子和那些原始變量相關(guān)，并分析該因子的隱含意義") print(paste(" 依據(jù)相對(duì)系數(shù)的絕對(duì)值大于在",loadings.abs.std,"原則，我們建議如下：")) tmp.x<-0#用于記錄因子載荷大于loadings.abs.std的所有原始變量的序列號(hào) dev.new()#新窗口# 畫出每個(gè)因子對(duì)應(yīng)各個(gè)變量的柱狀圖 par(mfrow=c(1,factors))#把窗口分為：1行3列 for(i in 1:factors){ tmp.x.factor<-0#用于記錄某一因子中因子載荷大于loadings.abs.std的原始變量的序列號(hào) print(paste(" 因子",i,"可以代表原始變量:")) print.con<-"" for(j in 1:length(names(data_frame))){ if(abs(loadings(sol)[,i][j])>loadings.abs.std){#loadings.abs.std為相對(duì)系數(shù)的絕對(duì)值的最低標(biāo)準(zhǔn) print(paste(" ",print.con,names(loadings(sol)[,i][j]),loadings(sol)[,i][j])) tmp.x<-c(tmp.x,j) tmp.x.factor<-c(tmp.x.factor,j) } } data1<-sol$loadings[,i] data1[-tmp.x.factor]<-0 data2<-sol$loadings[,i] data2[tmp.x.factor]<-0 barplot(data1,horiz=TRUE,main=paste("因子",i,"的載荷"),col="red",xlim=c(-1,1)) barplot(data2,horiz=TRUE,add=TRUE) } tmp.x<-tmp.x[-1] print(" 沒有被代表的原始變量有:") for(i in names(data_frame)[-as.numeric(names(table(tmp.x)))]){ print(paste(" ",i)) } tmp.x.table<-table(tmp.x) for(i in 1:length(tmp.x.table)){ if(tmp.x.table[i]>1){ print(paste(" Warings:原始變量",names(data_frame)[as.numeric(names(tmp.x.table[i]))],"被",tmp.x.table[i],"個(gè)因子共同代表了")) } } print("2 特殊值：") print("======================================================================") ssloadings<-sol$loadings[1,] for(i in 1:factors){ ssloadings[i]<-sum((sol$loadings[,i])^2) } var.sum<-length(names(data))#是每組xi數(shù)據(jù)標(biāo)準(zhǔn)化后(方差=1)的和=1*原始變量的個(gè)數(shù) tmp<-0 tmp.vector<-sol$loadings[1,]#每個(gè)因子對(duì)應(yīng)的累計(jì)貢獻(xiàn)比例 for(i in 1:factors){ print(paste(" 因子",i,"可以解釋所有原始變量X", round(10000*ssloadings[i]/var.sum)/100,"%的方差")) tmp<-(ssloadings[i]/var.sum)+tmp tmp.vector[i]<-tmp print(paste(" 因子1至",i,"累計(jì)可以解釋所有原始變量X", round(10000*tmp)/100,"%的方差")) print("") } print(" 請查看所有因子的累計(jì)方差貢獻(xiàn)比例，一般來說要大于80%，否則說明因子數(shù)目不足") dev.new()#新窗口#畫出累計(jì)方差貢獻(xiàn)比例 #barplot(tmp.vector,ylim=c(0,1),main="各個(gè)因子對(duì)整體方差的累計(jì)貢獻(xiàn)率ssloadings（特征值）") barplot(ssloadings/var.sum,ylim=c(0,1),main="各個(gè)因子對(duì)整體方差的貢獻(xiàn)率") lines(tmp.vector,col="blue",lwd=2) points(c(1:factors),tmp.vector) text(c(1:factors),tmp.vector+0.05,labels=round(tmp.vector*10000)/10000) abline(h=0.8,col="red") print("======================================================================") print("3 因子得分：") print(" 使用新產(chǎn)生的因子來表示原來的樣本") print(" 注意：每組因子對(duì)應(yīng)的樣本數(shù)據(jù)(即：每一列)已經(jīng)經(jīng)過了標(biāo)準(zhǔn)化：均值約為0，標(biāo)準(zhǔn)差約為1") print(sol$score) sol}##############系統(tǒng)聚類法：對(duì)新樣本進(jìn)行距離############################################hc<-function(data_frame,k.num){ d<-dist(scale(data_frame))#scale是標(biāo)準(zhǔn)化公式 hc<-hclust(d) dev.new() plclust(hc,hang=-1) re<-rect.hclust(hc,k=k.num,border="red")#劃分5個(gè)聚類，re[[i]]是第i個(gè)聚類包含的樣本id向量。 dev.new() hc.point<-data_frame[1:k.num,]#用于存儲(chǔ)每個(gè)聚類里個(gè)變量的平均值 for(i in 1:k.num){ for(j in 1:length(names(data_frame))){ hc.point[i,j]<-mean(data_frame[re[[i]],j]) } } #stars(hc.point+1+round(abs(min(hc.point))))#由于hc.point被標(biāo)準(zhǔn)化，所有有負(fù)數(shù)，無法使用星圖表示，現(xiàn)全體加一個(gè)數(shù)字，使不再有負(fù)數(shù)。 #dev.new() stars(hc.point+1+round(abs(min(hc.point))),full=F,draw.segments=T,key.loc=c(5,0.5),mar=c(2,0,0,0)) dev.new() stars(hc.point,full=F,draw.segments=T,key.loc=c(5,0.5),mar=c(2,0,0,0),main="不平移")#不知到底使用要平移 hc.point} |
測試程序：
補(bǔ)充數(shù)據(jù)applicant.csv
聯(lián)系客服
