#*************因子分析-R语言实现，函数库文件**************##****作者：oldlee11***************************************##****版本：v0.1*******************************************##****时间：2013-1-17*************************************##功能目标：原始数据变量x1,x2,x3....xn(全体记为X)。通过样本可以知道某些变量之间有相关性。# 则计算出新变量/因子f1,f2,f3....fm（m<n）（全体记为F），这新变量/因子F可以最大程度的表达原变量X# 由于新变量的个数m小于原始变量个数n，即降维了。#原理：AF+e=X# A叫因子载荷（loading）。意义：fi(某1个因子)和xi(某一个原变量)的相关系数，接近1表示fi和xi相关性强:aij=cov(xi,fj)# e叫特殊因子#其它术语变量：# 公因子方差：F（所有因子）解释xi（某一个原始变量）的方差百分比（贡献）# 特征值：fi（某一个因子）解释X（所有原始变量）的方差百分比（贡献）# 因子得分：在计算得出了A后，计算F内的样本数据。# 旋转：对因子载荷进行旋转，之后因子载荷各项大的越来越大，小的越来越小，便于划分因子和原始变量的关系。#****函数：factor()#****概要：因子分析法#****输入：# 名称 | 数据格式# data_frame | 欲分析的数据 ，数据框格式，最好带名称，不可以有factor分类数据# factors | 欲产生的因子个数小于原始数据的变量个数# scores | 是否进行因子得分# rotation | 是否进行旋转# loadings.abs.std | 用于从因子载荷中挑出那些原始变量和那些因子相关的标准（相对系数的绝对值的最低标准）#****输出：# factanal(..)产生的结果# factanal(..)$score并非数据框格式，需要as.data.frame(factanal(..)$score)转化一下才是数据框格式。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)#使用最大似然法进行的因子分析。 print("==========================") print("==== 因子分析结果如下 ====") print("==========================") print("模型：X=AF+e") print("======================================================================") print("1 A即因子载荷loadings:") print(" 每个数据(aij)表示了原始变量xi和因子变量fj的相关系数cov.值约接近+1或-1，约相关") for(i in 1:factors){ print(paste(" 因子",i,"同原始变量的相关性系数")) print(rev(sort(loadings(sol)[,i]))) } #### 给出因子和原始变量的可能关系#### print(" 您可以通过以上数据查看接近+1或者-1的数据，以说明某一因子和那些原始变量相关，并分析该因子的隐含意义") print(paste(" 依据相对系数的绝对值大于在",loadings.abs.std,"原则，我们建议如下：")) tmp.x<-0#用于记录因子载荷大于loadings.abs.std的所有原始变量的序列号 dev.new()#新窗口# 画出每个因子对应各个变量的柱状图 par(mfrow=c(1,factors))#把窗口分为：1行3列 for(i in 1:factors){ tmp.x.factor<-0#用于记录某一因子中因子载荷大于loadings.abs.std的原始变量的序列号 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为相对系数的绝对值的最低标准 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],"个因子共同代表了")) } } 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数据标准化后(方差=1)的和=1*原始变量的个数 tmp<-0 tmp.vector<-sol$loadings[1,]#每个因子对应的累计贡献比例 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,"累计可以解释所有原始变量X", round(10000*tmp)/100,"%的方差")) print("") } print(" 请查看所有因子的累计方差贡献比例，一般来说要大于80%，否则说明因子数目不足") dev.new()#新窗口#画出累计方差贡献比例 #barplot(tmp.vector,ylim=c(0,1),main="各个因子对整体方差的累计贡献率ssloadings（特征值）") barplot(ssloadings/var.sum,ylim=c(0,1),main="各个因子对整体方差的贡献率") 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(" 使用新产生的因子来表示原来的样本") print(" 注意：每组因子对应的样本数据(即：每一列)已经经过了标准化：均值约为0，标准差约为1") print(sol$score) sol}##############系统聚类法：对新样本进行距离############################################hc<-function(data_frame,k.num){ d<-dist(scale(data_frame))#scale是标准化公式 hc<-hclust(d) dev.new() plclust(hc,hang=-1) re<-rect.hclust(hc,k=k.num,border="red")#划分5个聚类，re[[i]]是第i个聚类包含的样本id向量。 dev.new() hc.point<-data_frame[1:k.num,]#用于存储每个聚类里个变量的平均值 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被标准化，所有有负数，无法使用星图表示，现全体加一个数字，使不再有负数。 #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} |
测试程序：
补充数据applicant.csv
联系客服
