因子分析(因子评价)

因子分析
一．因子分析原理
因子分析是根据相关性大小把原始变量进行分组，使得同组内的变量之间相关性高，而不同组的变量之间的相关性低。

每组变量代表一个基本结构（即公共因子），并用一个不可观测的综合变量来表示。

对于所研究的某一具体问题，原始变量分解为两部分之和。

一部分是少数几个不可观测的公共因子的线性函数，另一部分是与公共因子无关的特殊因子。

从全部计算过程来看作R 型因子分析与作Q 型因子分析都是一样的，只不过出发点不同，R 型从相关系数矩阵出发，Q 型从相似系数阵出发都是对同一批观测数据，可以根据其所要求的目的决定用哪一类型的因子分析
因子模型的性质：模型不受变量量纲的影响；因子载荷不是唯一的。

二．因子分析的数学模型
设有p 个指标，则因子分析数学模型为：
11111221221122221122
p p p p
p p p pp p X r Y r Y r Y X r Y r Y r Y X r Y r Y r Y
=+++⎧⎪
=+++⎪⎨
⎪
⎪=+++⎩ 其中，12,,
,p X X X 是已标准化的可观测的评价指标。

12,,,k F F F 出现在每个指标i X 的
表达式中，称为公共因子，公共因子是不可观测的，其含义要根据具体问题来解释。

i ε是
各个对应指标i X 所特有的因子，故称为特殊因子，它与公共因子之间彼此独立。

ij r 是指标i X 在公共因子j F 上的系数，称为因子载荷，因子载荷ij r 的统计含义是指标i X 在公共因子j F 上的相关系数，表示i X 与j F 线性相关程度。

用矩阵形式表示为：
X AF ε=+
其中12(,,
,)p X X X X '=，12(,,,)k F F F F '=，12(,,,)p εεεε'=，
1112
12122
212
m m p p pm r r r r r r A r
r r ⎛⎫
⎪
⎪
= ⎪ ⎪ ⎪⎝⎭
，A 称为因子载荷矩阵。

其统计含义是：
A 中的第i 行元素12,,,i i im r r r 说明了指标i X 依赖于各个公共因子的程度。

A 中第j 列元素12,,
,j j mj r r r 说明了公共因子j F 与各个指标的联系程度。

故常根据该
列绝对值较大的因子载荷所对应的指标来解释这个公共因子的实际意义。

A 中的第i 行元素12,,
,i i im r r r 的平方和2
21
m
i ij j h r ==∑称为指标i X 的共同度。

A 中第j 列元素12,,
,j j mj r r r 的平方和2
21
p
j ij i g r ==∑表示公共因子j F 对原始指标所提
供的方差贡献的总和，衡量各个公共因子的相对重要性。

称2
11p j
j ij i g r p p α===∑为公共因子j F 的方差贡献率，j α越大，公共因子j F 越重要。

三．因子分析的步骤
3.3.1
将原始变量数据进行标准化处理i Z =；
3.2 计算标准化指标的相关系数矩阵R ；
3.3 求解相关系数矩阵R 的特征向量()ij p p u u ⨯=和特征值120p λλλ≥≥
≥≥；
3.4 确定公共因子的个数，设为m 个，即选择特征值≥1的个数m 或根据累积方差贡献率
≥85％的准则所确定的个数m 为公共因子个数； 3.5
求解初始因子载荷矩阵()(ij p p p p A a u ⨯⨯==；
常用的方法有：主成分法、主轴因子法、极大似然法等。

本文用主成分法寻找公因子的方法如下：
设从相关矩阵出发求解主成分，设有p 个变量，则可以找出p 个主成分，将所得的p 个主成分由大到小排列，记为12,,
,p Y Y Y ，则主成分与原始变量之间有
11111221221122221122
p p
p p
p p p pp p Y r X r X r X Y r X r X r X Y r X r X r X
=+++⎧⎪
=+++⎪⎨
⎪
⎪=+++⎩ 其中ij r 是随机变量X 的相关矩阵的特征值所对应的特征向量的分量，特征向量之间正交，从X 到Y 的转换关系的可逆得到由Y 到X 的转换关系
11111221221122221122
p p
p p
p p p pp p X r Y r Y r Y X r Y r Y r Y X r Y r Y r Y
=+++⎧⎪
=+++⎪⎨
⎪
⎪=+++⎩ 只保留前m 个主成分，而把后面的p m -个主成分用特殊因子i ε代替，即
111112211221122222
1122m m m m p p p pm m p X r Y r Y r Y X r Y r Y r Y X r Y r Y r Y εεε=++++⎧⎪=++++⎪⎨
⎪
⎪=++++⎩
为了把i Y 转化为合适的公因子，需要把主成分i Y 变为方差为1的变量，故 令
i F =
ij a γ=则
111112211221122222
1122m m m m p p p pm m p X a F a F a F X a F a F a F X a F a F a F εεε=++++⎧⎪=++++⎪⎨
⎪
⎪=++++⎩
设样本相关系数矩阵R 的特征值为120p λλλ≥≥≥≥，其相应的标准正交特征向量
为12,,
,p γγγ，设m p <，则因子载荷矩阵A 的一个估计值为：
1211
2222ˆ(,m m m m
m p p pm u u u u u u A u u u λλλλγγγλλ⎛⎫
⎪
== ⎪
⎝ 共同度的估计为：22212ˆˆˆˆi
i i im h a a a
=+++。

3.6 建立因子模型
1
k
j ij j i i j Z a F a ε==+∑，1,2,
,i p =
其中12,,,k F F F 为公共因子，12(,,,)p εεεε=为特殊因子。

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

3.8 对初始因子载荷矩阵进行旋转 由于因子载荷矩阵不唯一，旋转变换可以是使初始因子载荷矩阵的每列或每行的元素的平方值趋于0或1，从而使得因子载荷矩阵结构简化，关系明确。

如果初始因子之间不相关，公共因子j F 的解释能力能够用其因子载荷平方的方差来度量时，则可采用方差极大正交旋转法；如果初始因子之间相关，则需要进行斜交旋转，通过旋转后，得到比较理想的新的因子载荷矩阵1()ij p k A r ⨯'=。

3.9 将公共因子变为变量的线性组合，得到因子得分函数
11221
k
i ij ij i i ip p j F Z Z Z Z ββββ===++
+∑，1,2,
,i m =
系数11B A R β-''=，i F ，ji Z 均为标准化的原始变量和公共因子。

因子得分函数的估计值为
11
11111ˆp m mp p b b X F A R X b b X -⎛⎫⎛⎫ ⎪⎪'==
⎪⎪ ⎪⎪⎝⎭⎝⎭
其中1A 为因子载荷矩阵，R 为原始变量的相关矩阵，X 为原始变量向量。

3.10 求综合评价值，即总因子得分估计值为
1
ˆˆm
i i i Z F ω==∑
其中1
i
i m
j
j λωλ
==
∑时第i 个公共因子i F 的归一化权重。

即：
⨯∑∑（各因子得分各因子所对应的方差贡献率）综合得分＝
各因子的方差贡献率
3.11 根据总因子得分估计值ˆZ
就可以对每个被评价的对象进行排名，从而进行比较。

四．因子分析的评价
4.1 首先在进行因子分析时，必须消除原始变量数据量纲和数量级的影响，所以需要对原始变量数据作转换。

常选用标准化变换i Z =
有些参考文献中也有说这样的标准化处
理仍然存在有不合理的地方，但是在实际应用中，为了简便，常选用上式进行变换。

4.2 在做因子分析之前，需要对原始变量间作相关性分析。

因为并不是所有的变量数据都是可以做因子分析的。

4.3因子分析适宜针对大样本容量做综合分析，对于小样本容量所做的分析不够准确。

一般要求样本容量大于指标个数的两倍。

4.4 不能简单地将初始因子载荷矩阵认为是主成分系数矩阵(特征向量矩阵)，否则会造成偏差。

所建立的综合评价函数1122m m F Y Y Y ααα=++
只是给出了一个排名，只是定性说明
这个函数包含了原始变量信息量的程度，并没有给出一个百分比等定量的度量。

五．典型案例
随着市场竞争的日益激烈，公司在人才选择方面更加注重人才的综合素质，并结合职位特定选择专门人才。

在本文中选取一家集生产与销售于一体的大公司在人才招聘中数据，从综合素质以及招聘职位来选择优秀的员工。

“华威”公司是一家集生产、销售为一体的大型国际著名公司。

现公司计划录用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 and Bartle tt's Te s t
.784645.317
105.000
Kaiser-Meyer-Olkin Measure of Sampling Adequacy.Approx. Chi-Square df Sig.
Bartlett's Test of Sphericity
通过观察上面的计算结果，可以知道，KMO 值为0.784，在较好的范围内；而巴特利球体检验的sig 值为0.00，拒绝原假设，说明相关矩阵并非单位矩阵，变量的相关系数较为显著。

3、求解特征根及相应特征向量
⏹Spss选项:Analyze－Data Reduction－Factor
⏹用Extraction，选择提取共因子的方法（如果是主成分分析，则选Principal
Components），
⏹用Rotation，选择因子旋转方法（如果是主成分分析就选None），
⏹用Scores计算因子得分，再选择Save as variables（因子得分就会作为变量存在
数据中的附加列上）和计算因子得分的方法（比如Regression）；要想输出Component
缓，解释能力不强。

因此我们提取5个公因子。

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

因子载荷矩阵
Com ponent Matrix
a .445.615.381-.103.098.584-.051-.028.287.748.110.340-.519.696-.183.617-.186.562.378-.111.796-.357-.291-.189.004.863-.188-.181-.078-.178.433-.581.343.456-.055.889-.042-.224-.217.032.367.793.100.074-.089.864.066-.096-.171-.173.872-.098-.252-.218.140.909-.033-.141.082-.071.914.032-.088.206-.109.711-.118.564-.220-.096.647.603
.108
-.022.070
求职信形式
外貌
专业能力
讨人喜欢能力自信心洞察力诚实
推销能力经验
驾驶水平事业心理解能力潜在能力交际能力适应性
12345Component
5 components extracted.
a.
通过因子载荷矩阵可以看出因子的意义并不是十分明确，为了对因子进行解释与说明，进行因子旋转，选取方差最大因子旋转方法，并保留因子得分。

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

最后，通过上面的因子选注我们的评价指标可以通过五个主要的因子来表示，分别为基本素质，工作经验素质，外在能力，专业素质和外表。

接下来计算各因子得分，并按照要求筛选优秀的应试者。

5、计算因子得分
各因子得分的计算公式为：
适应性
洞察力
自信心讨人喜欢能力
专业能力
外貌
求职信形式
*
025
.0
216
.0
*
249
.0
* 14 .0
*
005
.0
*
128
.0
*
112
.0
1
-
+
+
+-
-
-
-
=
f
适应性
洞察力
自信心讨人喜欢能力
专业能力
外貌
求职信形式
*
082
.0
209
.0
*
007
.0
* 062 .0
*
071
.0
*
033
.1
*
085
.0
5
+
+
-
+-
-
+
=
f
分别计算各应试者的五个因子得分，按照相对方差贡献率进行加权，得到最终各应试者的综合评价。

将综合得分进行排序，前六名为我们重点考虑的对象，编号分别为：10，1，2，22，14，
通过对排在前六名的应试者各因子得分分析可知：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_ro w-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;。