因子模型: X=μ + A*F* + ε
其中F=[(f1,f2,…,fm)]^T为公共因子向量,[ε=(ε1,ε2,…,εp)]^T为特殊因子向量,A=[(aij)]^(p×m)为因子载荷矩阵。
I.参数估计
为了建立因子模型,需要要得到因子载荷矩阵A=[(aij)]^(p×m)和特殊方差矩阵D=diag(σ1^2,σ2^2,…,σp^2)这两个参数的估计。
常用的参数估计方法有如下三种:主成分法、主因子法和极大似然法。
接下来会分别介绍以上三种方法具体方法,和综合三种方法的一个简便写法。
例. 12项智力指标的因子分析
研究者收集了40名学生的12项智力指标,分别为常识(x1)、类同(x2)、计算(x3)、词汇(x4)、理解(x5)、数字广度(x6)、常填图(x7)、图片排列(x8)、积木(x9)、拼图(x10)、译码(x11)和迷津(x12)。将原始数据经过标准化处理后,计算其相关系数矩阵,结果列在下表中。取m=2,试进行因子分析
#输入相关矩阵的数值
x <- c(
1.000,
0.6904 ,1.000,
0.4115 ,0.4511, 1.000,
0.4580, 0.7068, 0.4018, 1.000,
0.5535, 0.6620, 0.4122, 0.7119, 1.000,
0.3923, 0.6317, 0.4520, 0.4583, 0.5299, 1.000,
0.1415, 0.3009, 0.2025, 0.2665, 0.2480, 0.1590, 1.000,
0.0077, 0.0344, 0.1855, 0.1065, 0.0003, 0.1100, 0.3595, 1.000,
0.2385, 0.3523, 0.3646, 0.3644, 0.3388, 0.3982, 0.5004, 0.3314, 1.000,
0.0333, 0.1726, 0.1311, 0.1757, 0.1998, 0.0342, 0.5758, 0.1420, 0.2808, 1.000,
0.0898, 0.3878, 0.2041, 0.3191, 0.3186, 0.2914, 0.2537, 0.2025, 0.3971, 0.1468, 1.000,
0.2215, 0.2427, 0.4124, 0.2169, 0.1459, 0.0985, 0.4222, 0.2156, 0.5016, 0.2286, 0.0776, 1.000)
names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11", "X12")
R<-matrix(0, nrow=12, ncol=12, dimnames=list(names, names))
#生成相关系数矩阵R
for (i in 1:12){
for (j in 1:i){
R[i,j]<-x[(i-1)*i/2+j]; R[j,i]<-R[i,j]
}
}
1.主成分法
(需要设置的参数是R,因子个数m,后面会讲到m如何选取)
下面给出主成分法的R程序(factor.analy1.R)
factor.analy1<-function(S, m){
p<-nrow(S); diag_S<-diag(S); sum_rank<-sum(diag_S)
rowname<-paste("X", 1:p, sep="")
colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
eig<-eigen(S)
for (i in 1:m)
A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
h<-diag(A%*%t(A))
rowname<-c("SS loadings", "Proportion Var", "Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Principal Component Method")
list(method=method, loadings=A,
var=cbind(common=h, spcific=diag_S-h), B=B)
}
函数输入值S是样本方差阵或相关矩阵,m是主因子的个数,函数的输出值是列表形式,其内容有估计参数的办法(主成分法),因子载荷(loadings),共性方差和特殊方差,以及因子F对变量X的贡献、贡献率和累积贡献率。
#调用因子分析主成分法的函数
source("factor.analy1.R")
#显示结果.估计参数的方法为主成分法,loadings-因子载荷,var-共性方差和特殊方差,以及B-因子F对变量X的贡献、贡献率和累积贡献率
fa1<-factor.analy1(R, m=2); fa1
#协方差阵S的近似公式,误差平方和Q(m) (近似公式为E=S-A*A^T-D)
E1 <- R-fa1$loadings %*% t(fa1$loadings)-diag(fa1$var[,2])
sum(E1^2)
因子个数m的选取
#求特征值,对其求和
eigen(cor(R))
sum(eigen(cor(R))$values)
#选取满足 m个λ累加/所有λ累加 >= P0 的最小m,P0一般取[0.7,1)
(5.561644e+00 + 1.676901e+00 + 1.434965e+00) / sum(eigen(cor(R))$values)
#可以取m=3
#下面检验是否此时Q(m)最小
fa11 <- factor.analy1(R, m=3); fa11
#协方差阵S的近似公式,误差平方和Q(m) (近似公式为E=S-A*A^T-D)
E11 <- R-fa11$loadings %*% t(fa11$loadings)-diag(fa11$var[,2])
sum(E11^2)
结果看到,sum(E1^2)=1.060286 > sum(E11^2)=0.9550174。说明公因子个数m选择适当时,近似公式S的误差平方和Q(m)更优
2.主因子法
(需要设置的参数是R,因子个数m,特殊方差的估计值d,m值选取参考主成分法,d选取方法后面会讲到)
按照主因子法的思想编写相应的R程序:(factor.analy2.R)
factor.analy2<-function(R, m, d){
p<-nrow(R); diag_R<-diag(R); sum_rank<-sum(diag_R)
rowname<-paste("X", 1:p, sep="")
colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
kmax=20; k<-1; h <- diag_R-d
repeat{
diag(R)<- h; h1<-h; eig<-eigen(R)
for (i in 1:m)
A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
h<-diag(A %*% t(A))
if ((sqrt(sum((h-h1)^2))<1e-4)|k==kmax) break
k<-k+1
}
rowname<-c("SS loadings", "Proportion Var", "Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Principal Factor Method")
list(method=method, loadings=A,
var=cbind(common=h, spcific=diag_R-h), B=B, iterative=k)
}
函数输入值R是样本方差阵或相关矩阵,m是主因子的个数,d是特殊方差的估计值,函数的输出值是列表形式,其内容有估计参数的办法(主因子法),因子载荷(loadings),共性方差和特殊方差,以及因子F对变量X的贡献、贡献率和累积贡献率,以及求解的迭代次数。
相同数据,相关系数矩阵R,取公因子个m=2,特殊方差的估计值为0
#输入特殊方差var$spcific估计值,可以全部取0,下面会介绍怎么取合适的特殊方差估计值
d<-c(0,0,0,0,0,0,0,0,0,0,0,0)
#调用调用因子分析主因子法的函数
source("factor.analy2.R")
#显示结果.估计参数的方法为主成分法,loadings-因子载荷,var-共性方差和特殊方差,以及B-因子F对变量X的贡献、贡献率和累积贡献率,iterative-迭代次数
fa2<-factor.analy2(R, m=3, d); fa2
#近似公式S的误差平方和Q(m)
E2<- R-fa2$loadings %*% t(fa2$loadings)-diag(fa2$var[,2])
sum(E2^2)
用了13次迭代得到稳定解,再计算Q(m)
sum(E2^2)=0.3141111,优于主成分法
特殊方差估计值σi^2的常用选取方法
## 1.σi^2 = 1/rii,其中rii为R的逆矩阵的第i个对角线元素,此时Q(m)=sum(E21^2)
#R的逆矩阵R^-1
solve(R)
#取其对角线值,再求倒数
1 / diag(solve(R))
#将刚才的结果作为特殊方差估计值,我们来验证是否Q(m)会更优
d1 <- c(0.4113202,0.2159605,0.5974511,0.3610979,0.3659987,0.4522035,0.4673815,0.7639169,0.4743578,0.6385381,0.6627739,0.5743706)
#调用调用因子分析主因子法的函数
source("factor.analy2.R")
#显示结果.估计参数的方法为主成分法,loadings-因子载荷,var-共性方差和特殊方差,以及B-因子F对变量X的贡献、贡献率和累积贡献率,iterative-迭代次数
fa21 <- factor.analy2(R, m=3, d1); fa21
#近似公式S的误差平方和Q(m)
E21 <- R-fa21$loadings %*% t(fa21$loadings)-diag(fa21$var[,2])
sum(E21^2)
## 2.σi^2 = 1-hi^2,其中hi^2=max(j/i) |rij|
## 3.σi^2 = 1-hi^2,其中hi^2=1,此时σi^2全取0,此时Q(m)=sum(E2^2)
#### 这里R为相关矩阵,对角线元素全为1,其余元素都为0-1间的小数,所以方法2.和3.在这里是一样的
sum(E21^2) = 0.3106186 < sum(E2^2)=0.3141111,证明特殊方差估计值的选取方法,1.要优于2.、3.
3.极大似然法
(需要设置的参数是R,因子个数m,特殊方差的估计值d,m值选取参考主成分法,d值选取参考主因子法)
按照极大似然法的思想编写相应的R程序:(factor.analy3.R)
factor.analy3<-function(S, m, d){
p<-nrow(S); diag_S<-diag(S); sum_rank<-sum(diag_S)
rowname<-paste("X", 1:p, sep="")
colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
kmax=20; k<-1
repeat{
d1<-d; d2<-1/sqrt(d); eig<-eigen(S * (d2 %o% d2))
for (i in 1:m)
A[,i]<-sqrt(eig$values[i]-1)*eig$vectors[,i]
A<-diag(sqrt(d)) %*% A
d<-diag(S-A%*%t(A))
if ((sqrt(sum((d-d1)^2))<1e-4)|k==kmax) break
k<-k+1
}
rowname<-c("SS loadings","Proportion Var","Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Maximum Likelihood Method")
list(method=method, loadings=A,
var=cbind(common=diag_S-d, spcific=d),B=B,iterative=k)
}
函数输入值R是样本方差阵或相关矩阵,m是主因子的个数,d是特殊方差的估计值,函数的输出值是列表形式,其内容有估计参数的办法(主因子法),因子载荷(loadings),共性方差和特殊方差,以及因子F对变量X的贡献、贡献率和累积贡献率,以及求解的迭代次数。
相同数据,相关系数矩阵R,取公因子个m=2,特殊方差的估计值为:
#输入特殊方差var$spcific估计值(用上例中方法1.的结果d1)
d1 <- c(0.4113202,0.2159605,0.5974511,0.3610979,0.3659987,0.4522035,0.4673815,0.7639169,0.4743578,0.6385381,0.6627739,0.5743706)
#调用调用因子分析极大似然法的函数
source("factor.analy3.R")
#显示结果.估计参数的方法为主成分法,loadings-因子载荷,var-共性方差和特殊方差,以及B-因子F对变量X的贡献、贡献率和累积贡献率,iterative-迭代次数
fa3 <- factor.analy3(R, m=3, d1); fa3
#近似公式S的误差平方和Q(m)
E3 <- R-fa3$loadings %*% t(fa3$loadings)-diag(fa3$var[,2])
sum(E3^2)
sum(E3^2) = 0.3412492
4.综合以上三种方法
(method=“xxx”)
将上述3种方法结合在一起,并考虑主成分估计中介绍的因子个数m的选取方法,和在主因子法中介绍的特殊方差初始估计方法,编写相应的R程序
factor.analy.R
用一条函数,通过改变参数method=“xxx” , 可以更方便对比三种方法的结果
factor.analy<-function(S, m=0,
d=1/diag(solve(S)), method="likelihood"){
if (m==0){
p<-nrow(S); eig<-eigen(S)
sum_eig<-sum(diag(S))
for (i in 1:p){
if (sum(eig$values[1:i])/sum_eig>0.70){
m<-i; break
}
}
}
source("factor.analy1.R")
source("factor.analy2.R")
source("factor.analy3.R")
switch(method,
princomp=factor.analy1(S, m), #method=“princomp”时输入S,m=i两个参数
factor=factor.analy2(S, m, d), #method=“factor”时输入S,m=i,d=c(x,..,x)三个参数
likelihood=factor.analy3(S, m, d) #method=“likehood”时输入S,m=i,d=c(x,..,x)三个参数
)
}
函数输入样本方差矩阵S或样本相关矩阵R。因子个数m(缺省值由贡献率计算出m值)。特殊方差的初始估计d(缺省值为^σi方 = 1/rii)
计算因子载荷的方法,method=princomp采用主成分法,method=factor采用主因子法,method=likelihood(缺省值)采用极大似然法
函数输出就是采用前面介绍的三种方法的输出格式。
#使用factor.analy.R的实例:
source("factor.analy.R")
fa4 <- factor.analy(S=R,m=3,method = "princomp") ; fa4
#近似公式S的误差平方和Q(m)
E4 <- R-fa4$loadings %*% t(fa4$loadings)-diag(fa4$var[,2])
sum(E4^2) #可以看到,这里E4算出的Q(m)与E11算出的Q(m)是相同的
II.方差最大的正交旋转
某医院为了合理评价该院各月的医疗工作质量,搜集了3年有关门诊人次、出院人数、病床利用率、病床周转次数、平均住院天数、治愈好转率、病死率、诊断符合率、抢救成功率等9个指标数据,试采用因子分析法,探讨其综合评价指标体系。
使方差最大的因子载荷矩阵
先用三种方法之一计算的因子载荷估计矩阵,再用varimax()函数得到方差最大的因子载荷矩阵
#导入原始数据
hospital <- read.csv("hospital.csv",header=T)
#生成hospital表格的相关系数矩阵R
R <- cor(hospital)
for (i in 1:9){
for (j in 1:i){
R[i,j]<-x[(i-1)*i/2+j]; R[j,i]<-R[i,j]
}
}
#调用因子分析的特殊方差初始估计方法
source("factor.analy.R")
#以princomp方法为例
fa<-factor.analy(R, m=2, method="princomp")
vm1<-varimax(fa$loadings, normalize = F); vm1
因子分析的计算函数
事实上,在R软件中,提供了作因子分析计算的函数–factanal()函数,它可以从样本数据、样本的方差矩阵和相关矩阵出发对数据作因子分析,并可直接给出方差最大的载荷因子矩阵。
#显示factanal()函数的帮助页面,参数设置问题
?factanal()
#取公因子个数m=2,选用II中例子里的相关系数矩阵R,利用factanal函数得到fa结果
fa <- factanal(factors = 4,covmat = R)
#或者不用相关系数矩阵R,直接用csv格式文件:fa <- factanal(X=~.,factors=2,data=hospital)
#显示结果
fa
在上述信息中,call表示调用函数的方法,uniquenesses是特殊方差,loadings是因子载荷矩阵,其中Factor1,Factor2是因子,X1,X2,…,X9是对应的变量,SS loadings是公共因子对变量X的总方差贡献,Proportion Var是方差贡献率,Cumulative Var是累积方差贡献率。
IV.因子得分
回归法和加权最小二乘法
##导入原始数据
hospital <- read.csv("hospital.csv",header=T)
#相关矩阵特征值
eigen(cor(R))$values
sum(eigen(cor(R))$values[1:3])/sum(eigen(cor(R))$values)
#前3个因子的累积贡献率达到0.8134434,接下来选取因子个数为3
#不同方法计算因子得分
fa_1<-factanal(~., factors=3, data=hospital, scores="Bartlett") #加权最小二乘法
fa_2<-factanal(~., factors=3, data=hospital, scores="regression") #回归法
fa_1;fa_2
#画出各组在第a、第b公共因子下的散点图
plot(fa$scores[, 1:2], type="n"); text(fa$scores[,1], fa$scores[,2]) #第一、第二公共因子下的散点图
plot(fa$scores[, c(1,3)], type="n"); text(fa$scores[,1], fa$scores[,3]) #第一、第三公共因子下的散点图
上面是采用回归法,也可以使用加权最小二乘法来画图。由散点图,可以直观选出偏向哪个公共因子的组别。
根据选项factors=4的设定,3个潜在因子被保留,前3个因子的累积贡献率达到81.3%,上式为全部变量在3个潜在因子F1-F3上的因子载荷矩阵
例如:x1由4个因子表达的式子为:
x1=0.447*F1 + 0.519*F2 + -0.101*F4
从矩阵上看,因子1在多数原始指标上均有较大的载荷,因子2在x1(门诊人次)、x2(出院人数)、x3(病床利用率)、x4(病床周转次数)上有较大的载荷,因子3在x2(出院人数)、x5(平均住院天数)、x6(治愈好转率)上有较大的载荷。除因子1可以认定为综合因子外,其他3个因子意义不明显。
因子旋转
fa_1<-factanal(~., factors=3, data=hospital, scores="Bartlett") #加权最小二乘法
vm <- varimax(fa_1$loadings,normalize = F) ; vm
经过因子旋转处理,3个潜在因子在9个原始指标上的因子载荷矩阵如上表所示。
对该因子载荷进行分析,可看出:因子F1在x1(门诊人次)、x2(出院人数)、x5(平均住院天数)、x8(诊断符合率)、x9(抢救成功率)上因子载荷较大;F2在x3(病床利用率)、x4(病床周转次数)上的因子载荷较大;F3在x6(治愈好转率)、x7(病死率)上的因子载荷较大
我们可以推出:因子F1反映了该医院医疗工作质量各方面的情况,为综合因子;F2反映了病床利用情况;F3反映了医疗水平的高低
将旋转后的因子载荷与主成分分析的因子载荷矩阵比较可知:因子旋转后,除F1的因子载荷仍分布多数指标上外,其他2个因子的载荷明显地集中到少数指标上,说明旋转对因子载荷起到明显的分离作用,使得各因子解释的变量更加清晰。
数据分析咨询请扫描二维码
CDA数据分析师在中国航信高科技产业园进行了面向测试度量的数据分析培训课程,培训人数近2 ...
2024-05-01CDA数据分析师走进深圳迈瑞生物医疗电子股份有限公司,在迈瑞总部展开了为期两天的培训,本次课程参训人员线上及线下近百人, ...
2024-05-01CDA数据分析师在合肥市对合肥阳光新能源科技有限公司开展了为期8天的企业内训。 合肥阳光新能源科技 ...
2024-05-01CDA数据分析师走进海尔大学,进行了《数据治理与数据中台建设的道与术》专题培训,培训现场爆满,近百人参加了此次培训。 ...
2024-05-01在中国银行苏州分行培训中心开始数据分析师培训,此次培训课程共10天内容,包括Excel、MySQL、概率论与数理统计、SPSS等内容, ...
2024-05-01从实际的业务需求出发,结合行业的典型应用特点,围绕实际的商业问题,探讨数据挖掘、机器学习模型在金融领域的应用,包括获客、信用评分、细分画像、交叉销售、反欺诈、违规识别、时序预测、运筹优化、流程挖掘九个方面,形成 ...
2024-05-01本次培训课程为线上+线下的模式,由于学员编程能力不一、部分学员没有编程基础,故提供统计学、python基 ...
2024-05-01华夏银行信用卡中心-机器学习培训 1、课程亮点 取材于业界一流企业和顶级咨询公司的行业实践;已经被证明是人人 ...
2024-05-01主 题:数据中台建设及数据分析应用主题分享 1. 数据中台市场洞察 2. 主流数据中台产品比较 3. 某企业数据中 ...
2024-05-01围绕“数据驱动”战略,全力打造我行 300 人数字化人才梯队,着力培养数字化管理人才、大数据专业团队 ...
2024-05-01在当今数据驱动的商业环境中,数据分析成为了企业决策的重要依据。通过对大量数据的收集、处理和分析,企业能够更好地理解市场 ...
2024-04-29在人工智能(AI)的世界里,提示词(Prompt)是一种强大的工具,它能够引导AI按照用户的需求产生特定的输出。本文将深入探讨AI ...
2024-04-29CDA立足未来职场,拓展前沿视野——对外经贸大学保险学院举办“三全育人大讲堂”分享行业最新动态。 ...
2024-04-294月2日,CDA数据分析师创始发起人兼协会理事长赵坚毅博士受邀在浙江万里学院举办了一场以“数字化能力在职场中的作用” ...
2024-04-29随机森林(Random Forests)现在机器学习中比较火的一个算法,是一种基于Bagging的集成学习方法,能够很好地处理分类和回归的问 ...
2022-12-23方差分析是数据分析中常用的一种统计分析方法,接下来让我们简单了解一下方差分析的基本思想和原理吧。 方差分析(Analysis ...
2022-12-23来源:关于数据分析与可视化 关于streamlit-aggrid 数据排序 表格样式的调整 数据 ...
2022-08-03作者:麦叔 定义 「把上面晦涩的概念汇成一句话就是:」 ❝ 回调函数就是一个被作为参 ...
2022-08-03现今,高学历人群日益增多,物以稀为贵的高学历光环淡去。无论本科生还是研究生,甚至博士生,求职竞争力都大不如前,就业压力越来越大。
2022-06-01某家企业10个人面试,有9个本科生……如何脱颖而出,除得体的举止和良好的沟通力外,证书成重要筹码,这也是很多人考证的关键所在。
2022-04-14