因子分析(因子评价)

因子分析

一．因子分析原理

    因子分析是根据相关性大小把原始变量进行分组，使得同组内的变量之间相关性高，而不同组的变量之间的相关性低。每组变量代表一个基本结构（即公共因子），并用一个不可观测的综合变量来表示。对于所研究的某一具体问题，原始变量分解为两部分之和。一部分是少数几个不可观测的公共因子的线性函数，另一部分是与公共因子无关的特殊因子。

从全部计算过程来看作R型因子分析与作Q型因子分析都是一样的，只不过出发点不同，R型从相关系数矩阵出发，Q型从相似系数阵出发都是对同一批观测数据，可以根据其所要求的目的决定用哪一类型的因子分析

因子模型的性质：模型不受变量量纲的影响；因子载荷不是唯一的。

二．因子分析的数学模型

设有个指标，则因子分析数学模型为：

其中，是已标准化的可观测的评价指标。出现在每个指标的表达式中，称为公共因子，公共因子是不可观测的，其含义要根据具体问题来解释。是各个对应指标所特有的因子，故称为特殊因子，它与公共因子之间彼此。是指标在公共因子上的系数，称为因子载荷，因子载荷的统计含义是指标在公共因子上的相关系数，表示与线性相关程度。

    用矩阵形式表示为：

其中，，，，称为因子载荷矩阵。

其统计含义是：

中的第行元素说明了指标依赖于各个公共因子的程度。

中第列元素说明了公共因子与各个指标的联系程度。故常根据该列绝对值较大的因子载荷所对应的指标来解释这个公共因子的实际意义。

中的第行元素的平方和称为指标的共同度。

中第列元素的平方和表示公共因子对原始指标所提供的方差贡献的总和，衡量各个公共因子的相对重要性。称为公共因子的方差贡献率，越大，公共因子越重要。

三．因子分析的步骤

3.3.1 将原始变量数据进行标准化处理；

3.2 计算标准化指标的相关系数矩阵；

3.3 求解相关系数矩阵的特征向量和特征值；

3.4 确定公共因子的个数，设为个，即选择特征值1的个数或根据累积方差贡献率85％的准则所确定的个数为公共因子个数；

3.5 求解初始因子载荷矩阵；

    常用的方法有：主成分法、主轴因子法、极大似然法等。本文用主成分法寻找公因子的方法如下：

    设从相关矩阵出发求解主成分，设有个变量，则可以找出个主成分，将所得的个主成分由大到小排列，记为，则主成分与原始变量之间有

其中是随机变量的相关矩阵的特征值所对应的特征向量的分量，特征向量之间正交，从到的转换关系的可逆得到由到的转换关系

只保留前个主成分，而把后面的个主成分用特殊因子代替，即

为了把转化为合适的公因子，需要把主成分变为方差为1的变量，故

令    ，

则                    

设样本相关系数矩阵的特征值为，其相应的标准正交特征向量为，设，则因子载荷矩阵的一个估计值为：

共同度的估计为：。

3.6 建立因子模型

，

其中为公共因子，为特殊因子。

3.7 对公共因子进行重新命名，并解释公共因子的实际含义

    当初始因子载荷矩阵难以对公共因子的实际意义作出解释时，先要对作方差极大正交旋转，然后再根据旋转后所得的正交因子载荷矩阵作出解释，即根据指标的因子载荷绝对值的大小，值的正负符号来说明公共因子的意义。

3.8 对初始因子载荷矩阵进行旋转

    由于因子载荷矩阵不唯一，旋转变换可以是使初始因子载荷矩阵的每列或每行的元素的平方值趋于0或1，从而使得因子载荷矩阵结构简化，关系明确。如果初始因子之间不相关，公共因子的解释能力能够用其因子载荷平方的方差来度量时，则可采用方差极大正交旋转法；如果初始因子之间相关，则需要进行斜交旋转，通过旋转后，得到比较理想的新的因子载荷矩阵。

3.9 将公共因子变为变量的线性组合，得到因子得分函数

，

系数，，均为标准化的原始变量和公共因子。因子得分函数的估计值为

其中为因子载荷矩阵，为原始变量的相关矩阵，为原始变量向量。

3.10 求综合评价值，即总因子得分估计值为

其中时第个公共因子的归一化权重。即：

3.11 根据总因子得分估计值就可以对每个被评价的对象进行排名，从而进行比较。

四．因子分析的评价

4.1 首先在进行因子分析时，必须消除原始变量数据量纲和数量级的影响，所以需要对原始变量数据作转换。常选用标准化变换。有些参考文献中也有说这样的标准化处理仍然存在有不合理的地方，但是在实际应用中，为了简便，常选用上式进行变换。

4.2     在做因子分析之前，需要对原始变量间作相关性分析。因为并不是所有的变量数据都是可以做因子分析的。

4.3因子分析适宜针对大样本容量做综合分析，对于小样本容量所做的分析不够准确。一般要求样本容量大于指标个数的两倍。

4.4 不能简单地将初始因子载荷矩阵认为是主成分系数矩阵(特征向量矩阵)，否则会造成偏差。所建立的综合评价函数只是给出了一个排名，只是定性说明这个函数包含了原始变量信息量的程度，并没有给出一个百分比等定量的度量。

五．典型案例

随着市场竞争的日益激烈，公司在人才选择方面更加注重人才的综合素质，并结合职位特定选择专门人才。在本文中选取一家集生产与销售于一体的大公司在人才招聘中数据，从综合素质以及招聘职位来选择优秀的员工。

“华威”公司是一家集生产、销售为一体的大型国际著名公司。现公司计划录用6名的员工。经过初选，公司对48位应聘者进行面试，面试共有15项指标，这15项指标分别是：求职信的形式(FL)、外貌(APP)、专业能力(AA)、讨人喜欢(LA)、自信心(SC)、洞察力(LC)、诚实(HON)、推销能力(SMS)、经验(EXP)、驾驶水平(DRV)、事业心(AMB)、理解能力(GSP)、潜在能力(POT)、交际能力(KJ)和适应性(SUIT)。每项指标的分数是从0分到10分，0分最低，10分最高。每位求职者的15项指标的得分在文件（应聘者得分记录.xls）中。试从综合素质选出6名优秀员工，若将这6名员工分别分配到管理、销售和生产部门各2名，指出合理的分配方案。

（二）、分析过程详解

1、数据标准化

由于数据均为在面试中的打分成绩，量纲相同，并且观察数据的分布，并无异常值的出现，因此数据没有必要进行标准化，可以直接进行分析。

2、建立相关系数矩阵

利用SPSS软件，correlate功能计算相关系数矩阵，计算皮尔森相关系数并进行卡方双尾检验，可以看出变量间存在这很大的相关性。

进行相关系数矩阵检验——KMO测度和巴特利特球体检验：

KMO值：0.9以上非常好；0.8以上好；0.7一般；0.6差；0.5很差；0.5以下不能接受。巴特利特球体检验原假设H0:相关矩阵为单位阵

通过观察上面的计算结果，可以知道，KMO值为0.784，在较好的范围内；而巴特利球体检验的sig值为0.00，拒绝原假设，说明相关矩阵并非单位矩阵，变量的相关系数较为显著。

3、求解特征根及相应特征向量

⏹Spss选项:Analyze－Data Reduction－Factor

⏹用Extraction，选择提取共因子的方法（如果是主成分分析，则选Principal Components），

⏹用Rotation，选择因子旋转方法（如果是主成分分析就选None），

⏹用Scores计算因子得分，再选择Save as variables（因子得分就会作为变量存在数据中的附加列上）和计算因子得分的方法（比如Regression）；要想输出Component Score Coefficient Matrix表，就要选择Display factor score coefficient matrix；

输出结果如下：

碎石图：

通过此图可以明显看出前五个因子可以解释大部分的方差，到第六个因子以后，线逐渐平缓，解释能力不强。因此我们提取5个公因子。

方差贡献率

选择5个公因子，从方差贡献率可以看出，其中第一个公因子解释了总体方差的50.092%，四个公因子的累计方差贡献率为86.42%,可以较好的解释总体方差。

因子载荷矩阵

通过因子载荷矩阵可以看出因子的意义并不是十分明确，为了对因子进行解释与说明，进行因子旋转，选取方差最大因子旋转方法，并保留因子得分。

4、因子旋转

旋转后的因子载荷矩阵：

通过上表中旋转的结果，我们可以看出第一个公因子在自信心，洞察力，推销能力，驾驶水平，事业心，理解能力，潜在能力上有较大的载荷，可以将其命名为基本素质；第二个因子在求职信形式，经验，适应性上有较大的载荷，可以理解为工作经验素质；第三个因子在讨人喜欢能力，诚实，交际能力上有较大的载荷，可以命名为外在能力；第四个因子在专业能力上载荷较大，但在交际能力上的载荷为负相关，也从侧面反映了专业能力较强的人在交际上有一定的欠缺，这和目前一部分高校毕业生书本专业知识较强，但日常待人接物能力较差的现象相吻合，将其命名为专业素质；第五个因子仅在外貌上有较大的载荷，可以将其命名为外表。

最后，通过上面的因子选注我们的评价指标可以通过五个主要的因子来表示，分别为基本素质，工作经验素质，外在能力，专业素质和外表。接下来计算各因子得分，并按照要求筛选优秀的应试者。

5、计算因子得分

各因子得分的计算公式为：

分别计算各应试者的五个因子得分，按照相对方差贡献率进行加权，得到最终各应试者的综合评价。

综合得分：

编号	综合得分	编号	综合得分	编号	综合得分
1	0.9925	17	0.2369	33	-0.0634
2	0.9315	18	0.1480	34	-0.5107
3	0.7912	19	0.2378	35	-0.4718
4	0.6150	20	0.2073	36	-0.3979
5	0.6968	21	-0.0711	37	-0.1255
6	0.7259	22	0.8704	38	-0.2785
7	0.6070	23	0.0993	39	-0.5476
8	0.6007	24	0.0325	40	-0.6171
9	0.4299	25	-0.2851	41	-0.9400
10	1.0461	26	0.2296	42	-0.78
11	0.4309	27	-0.3275	43	-1.1902
12	0.3623	28	-0.0763	44	-1.0797
13	0.5168	29	0.2112	45	-0.9281
14	0.8424	30	-0.3032	46	-0.9682
15	0.3403	31	-0.0685	47	-0.9455
16	0.2117	32	-0.37	48	-1.0527

将综合得分进行排序，前六名为我们重点考虑的对象，编号分别为：10，1，2，22，14，3。我们重点分析这六个人的各方面得分，从而指定较为合理的部门分配。

编号	基本素质		工作经验素质	外在能力	专业素质	外表	综合得分
10	2.04521		-0.50527	-1.61184	1.493955	-0.28452	1.046139
1	1.147623		1.552198	1.040729	0.4904	-1.42199	0.992542
2	1.009694		1.585842	1.165884	0.556039	-1.53665	0.931513
22	1.960285		-0.547	-2.70356	1.486739	-0.22207	0.870439
14	1.629787		-0.83238	-0.90829	1.651657	-0.36539	0.84243
3	0.833452		1.442163	0.168301	0.429408	0.363725	0.7911
6	1.121359	0.401301		0.698105	-0.61087	-0.170	0.725859

通过对排在前六名的应试者各因子得分分析可知：10号和22号在基本素质包括自信心，洞察力，推销能力，驾驶水平，事业心，理解能力，潜在能力上的得分最高，这也是他最大的优势所在，在专业素质上的得分也很高，说明此人专业水平和好，但是在外在能力也就是交际能力，讨人喜欢能力以及第五因子外表上得分均为负值，说明此人的外在条件并不是很好；另外此人的工作经验素质的得分也为负值，其中包括求职信形式，经验，适应性等，说明欠缺工作经验，此人很有可能为刚毕业的应届生，基本素质以及专业知识过硬，但是社会生活能力较差。在本次所招收的职位中，生产部门的研发职位要求专业能力强，对于经验及外在条件的要求较低，10号和22号比较适合安排在生产部门。

1号和2号除了外表为负值以外，其他因子得分均为正值，此人外在能力和工作经验素质的得分较高，说明此人有一定的工作经验，社交能力较强。本次招聘中，销售部门对于员工的基本要求就是工作经验以及较强的社交能力，因此1号比较适合安排在销售部门。

3号各方面的得分均为正值，说明比较全面，其中工作经验素质得分最高，说明此人工作经验丰富并且适应性较强；管理部门要求员工素质要全面均衡，工作经验比较重要，所以3号适合安排在管理部门。

14号的情况与10号和22号类似，由于我们现在需要管理部门的员工，从各项得分来说，14号不太适合；所以我们后移一位，选择6号安排进管理部门，虽然6号在外表和专业素质上的得分为负值，但这对于管理部门来说并不是最重要的。6号在因子1即基本素质的得分高于3号，自信心，洞察力，推销能力，事业心，理解能力，潜在能力等对于管理阶层才是最重要的。

（三）、主要结论

通过上面对应试者的分析我们可以看出，作为应届的学生，在找工作时，公司重点考察的方面，以及我们所欠缺的。如前述分析中大部分的应试者专业能力较强但是待人接物的能力明显较差，这也是应届生求职中的主要障碍之一；还有应届生缺乏工作经验，适应性较差，这也是困扰大学生求职的问题之一。作为应届生可以从中找到自身的弱点加以改善；同时还要注意自身的优点，在求职中结合职位要求有的放矢，才能实现较好的就业。

公司在选拔人才时不能笼统的只看评价打分的平均分或者总分，应该从公司的实际出发，从招聘职位的具体要求出发，选择最适合企业发展的人才。在多指标的评价体系中，因子分析可以将空间降维，实现以较少的指标来解释大部分自变量的信息，从而简化分析步骤，评价起来更清晰简便。

六．MATLAB程序代码

clear all;

DATA=load('D:0.m');

DATA=double(DATA);

DATA=DATA';

TESTDATA=load('D:14f.m');

TESTDATA=double(TESTDATA);

% DATA=load('D:正常.txt');

% DATA=double(DATA);

% DATA=DATA(:,3:12);

% TESTDATA=load('D:异常.txt');

% TESTDATA=double(TESTDATA);

% TESTDATA=TESTDATA(:,3:12);

[Kp,T2]=tztq(DATA,TESTDATA);

function [contribution,T2,SPE,t2cl,s_cl] =  PCA_model(Xtrain,Xtest)

X_mean = mean(Xtrain);                

X_std = std(Xtrain);                   

[X_row ,X_col]= size(Xtrain); 

for i = 1:X_col

    Xtrain(:,i) = (Xtrain(:,i)-X_mean(i))./X_std(i);

    Xtest(:,i) = (Xtest(:,i)-X_mean(i))./X_std(i);

end

[U,S,V]=svd(Xtrain./sqrt(size(Xtrain,1)-1),0);        

D= S^2;

lamda=diag(D);

num_pc=1;

while sum(lamda(1:num_pc))/sum(lamda)<0.9

    num_pc=num_pc+1;

end

D=diag(lamda);

P=V(:,1:num_pc);

[a,b]=size(Xtest);

[r,y]=size(P*P');

I=eye(r,y);

e=Xtest*(I-P*P');

for i=1:a  

T2(i)=Xtest(i,:)*P*inv(D(1:num_pc,1:num_pc))*P'*Xtest(i,:)';

end

for l=1:a

    SPE(l)=e(l,:)*e(l,:)';

end

 for j=1:b

    contribution(j)=(norm(e(:,j)))^2;

 end

 

 t2cl=num_pc*(X_row-1)*(X_row+1)*icdf('f',0.99,num_pc,X_row-num_pc)/(X_row*(X_row-num_pc));

for i=1:3

     theta(i)=trace((D(num_pc+1:X_col,num_pc+1:X_col))^i);        

end

% 另一种SPE控制线算法

% h=(theta(1)^2)/theta(2);

% g=theta(2)/theta(1);

% conf=0.95;                              

% df=round(h);                           

% delta2a1=g*pinv(df,2);

h0=1-2*theta(1)*theta(3)/(3*theta(2)^2);

ca=icdf('norm',0.99,0,1);

s_cl=theta(1)*(ca*sqrt(2*theta(2)*h0^2)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0);

clear all;

X1=load('D:0.m');

Xtrain=X1';

Xtrain=double(Xtrain);

X2=load('D:14f.m');

Xtest=X2;

Xtest=double(Xtest);

% X1=load('D:正常br.txt');

% Xtrain=X1(:,3:62);

% Xtrain=double(Xtrain);

% X2=load('D:异常br.txt');

% Xtest=X2(:,3:62);

% Xtest=double(Xtest);

[contribution,T2,SPE,t2cl,s_cl]=PCA_model(Xtrain,Xtest);

    figure

    x=size(Xtest,1);

    plot(1:x,SPE,'k');

    hold on;

    plot(1:x,s_cl,'r-');

    title('SPE');

    hold off;

    figure

    plot(1:x,T2,'K');

    title('T2');

    hold on;

    plot(1:x,t2cl,'r-');

    hold off;

    figure

    bar(contribution,'group')

    title('贡献图');

function [Kp,T2]=tztq(ax,ay)

[Nx]=size(ax);

mean_X =  mean(ax);

axb=ax;

std_X=std(ax);

ax=ax-mean_X(ones(Nx,1),:);

std_X(find(std_X==0))=1;%数据预处理

ax=ax./std_X(ones(Nx,1),:);

c=10000;

% gama=0.05;

% ni=1;

% F1=ax(1,:);

% F=F1';

% for ib=2:Nx

%     for i=1:ni

%         for j=1:ni

%           

%   batar1(ib).block(i,j)=exp(-norm(ax(i,:)-ax(j,:))^2/c);

%         end

%         batar2(i,1)=exp(-norm(ax(i,:)-ax(ib,:))^2/c);

%         batar3(1,i)=exp(-norm(ax(ib,:)-ax(i,:))^2/c);

%     end

%     s1=exp(-norm(ax(ib,:)-ax(ib,:))^2/c);

%     batar= batar3(1,i)*inv(batar1(ib).block)* batar2(i,1);

%     s=(s1- batar)/s1;

% if s>sqrt(gama)

%       ni=ni+1;

%       F=[F ax(ib,:)'];

%     end

%     

%   

% end

% AX=F'%训练集基的提取结束

  [N]=size(ax,1);

    for i=1:N

        for j=1:N       

  K(i,j)=exp(-norm(ax(i,:)-ax(j,:))^2/c);%求核矩阵   

        end

    end

    n1=ones(N,N);

    N1=1/N*n1;

    Kp=K-N1*K-K*N1+N1*K*N1;

    [u,s,v]=svd(Kp/N);

    lamda=s;

    P=v;

    lamda=diag(lamda);

    B=length(find(lamda>1e-10));   %求非零的特征值个数

    %求主元个数

npc=1;

while sum(lamda(1:npc))/sum(lamda(1:B))<0.9

    npc=npc+1;

end

npc

%求特征空间有效维数

DimFS=1;

while sum(lamda(1:DimFS))/sum(lamda(1:B))<=0.99

    DimFS=DimFS+1;

end

lamda=diag(lamda);

for i=1:B

    

%     P(:,i)=P(:,i)/norm(P(:,i)*s(i,i));

P(:,i)=P(:,i)/(norm(P(:,i))*sqrt(N*lamda(i,i)));

end

[Ny]=size(ay,1);

mean_X =mean(axb);

std_X = std(axb);

[num_sample] = Ny;    

ay = ay-mean_X(ones(num_sample,1),:);   

ay = ay./std_X(ones(num_sample,1),:); 

% mean_y =  mean(ay);

% std_y=std(ay);

% ay =  ay-mean_y(ones(Ny,1),:);

% std_y(find(std_y==0))=1;%数据处理

% ay =  ay./std_y(ones(Ny,1),:);

for i=1:Ny

   for j=1:N          

       Ky(i,j)=exp(-norm(ay(i,:)-ax(j,:))^2/c);

   end

end

    t1=ones(1,N);

    t11=1/N*t1;

    for i=1:Ny

        kp1(i,:)= Ky(i,:)-t11*K- Ky(i,:)*N1+t11*K*N1;

    end

for i=1:Ny

   for k=1:B

       t(i,k)=P(:,k)'*kp1(i,:)';

   end

end

  %    求T2,SPE

  %    covtyb=inv(t'*t);

    for i=1:Ny

        T2(i)=t(i,1:npc)*inv(lamda(1:npc,1:npc))*t(i,1:npc)';   %也可以

%       SPE(i)=t(i,1:npc)*t(i,1:npc)';

%       T2(1,i)=t(i,1:npc)*(covtyb(1:npc,1:npc))*t(i,1:npc)'; 

        SPE(i)=t(i,(npc+1):B)*t(i,(npc+1):B)';

    end 

    %T2,SPE控制线

    t2cl=npc*(N-1)*(N+1)*icdf('f',0.99,npc,N-npc)/(N*(N-npc));

    for i=1:3

     theta(i)=trace((lamda(npc+1:DimFS,npc+1:DimFS))^i);        

    end

    h0=1-2*theta(1)*theta(3)/(3*theta(2)^2);

    ca=icdf('norm',0.99,0,1);

    s_cl=theta(1)*(ca*sqrt(2*theta(2)*h0^2)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0);

     figure

     plot(1:Ny,T2,'k');

     hold on;

     plot(1:Ny,t2cl,'r');

     title('T2');

     hold off;

     figure

     plot(1:Ny,SPE,'k')

     hold on;

     plot(1:Ny, s_cl,'r');

     title('SPE');

     hold off;