
因子模型: 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如何選?。?br />
下面給出主成分法的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)更優(yōu)
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次迭代得到穩(wěn)定解,再計算Q(m)
sum(E2^2)=0.3141111,優(yōu)于主成分法
特殊方差估計值σ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)會更優(yōu)
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.要優(yōu)于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=“l(fā)ikehood”時輸入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.方差最大的正交旋轉
某醫(yī)院為了合理評價該院各月的醫(yī)療工作質量,搜集了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()函數,它可以從樣本數據、樣本的方差矩陣和相關矩陣出發(fā)對數據作因子分析,并可直接給出方差最大的載荷因子矩陣。
#顯示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反映了該醫(yī)院醫(yī)療工作質量各方面的情況,為綜合因子;F2反映了病床利用情況;F3反映了醫(yī)療水平的高低
將旋轉后的因子載荷與主成分分析的因子載荷矩陣比較可知:因子旋轉后,除F1的因子載荷仍分布多數指標上外,其他2個因子的載荷明顯地集中到少數指標上,說明旋轉對因子載荷起到明顯的分離作用,使得各因子解釋的變量更加清晰。
數據分析咨詢請掃描二維碼
若不方便掃碼,搜微信號:CDAshujufenxi
LSTM 模型輸入長度選擇技巧:提升序列建模效能的關鍵? 在循環(huán)神經網絡(RNN)家族中,長短期記憶網絡(LSTM)憑借其解決長序列 ...
2025-07-11CDA 數據分析師報考條件詳解與準備指南? ? 在數據驅動決策的時代浪潮下,CDA 數據分析師認證愈發(fā)受到矚目,成為眾多有志投身數 ...
2025-07-11數據透視表中兩列相乘合計的實用指南? 在數據分析的日常工作中,數據透視表憑借其強大的數據匯總和分析功能,成為了 Excel 用戶 ...
2025-07-11尊敬的考生: 您好! 我們誠摯通知您,CDA Level I和 Level II考試大綱將于 2025年7月25日 實施重大更新。 此次更新旨在確保認 ...
2025-07-10BI 大數據分析師:連接數據與業(yè)務的價值轉化者? ? 在大數據與商業(yè)智能(Business Intelligence,簡稱 BI)深度融合的時代,BI ...
2025-07-10SQL 在預測分析中的應用:從數據查詢到趨勢預判? ? 在數據驅動決策的時代,預測分析作為挖掘數據潛在價值的核心手段,正被廣泛 ...
2025-07-10數據查詢結束后:分析師的收尾工作與價值深化? ? 在數據分析的全流程中,“query end”(查詢結束)并非工作的終點,而是將數 ...
2025-07-10CDA 數據分析師考試:從報考到取證的全攻略? 在數字經濟蓬勃發(fā)展的今天,數據分析師已成為各行業(yè)爭搶的核心人才,而 CDA(Certi ...
2025-07-09【CDA干貨】單樣本趨勢性檢驗:捕捉數據背后的時間軌跡? 在數據分析的版圖中,單樣本趨勢性檢驗如同一位耐心的偵探,專注于從單 ...
2025-07-09year_month數據類型:時間維度的精準切片? ? 在數據的世界里,時間是最不可或缺的維度之一,而year_month數據類型就像一把精準 ...
2025-07-09CDA 備考干貨:Python 在數據分析中的核心應用與實戰(zhàn)技巧? ? 在 CDA 數據分析師認證考試中,Python 作為數據處理與分析的核心 ...
2025-07-08SPSS 中的 Mann-Kendall 檢驗:數據趨勢與突變分析的有力工具? ? ? 在數據分析的廣袤領域中,準確捕捉數據的趨勢變化以及識別 ...
2025-07-08備戰(zhàn) CDA 數據分析師考試:需要多久?如何規(guī)劃? CDA(Certified Data Analyst)數據分析師認證作為國內權威的數據分析能力認證 ...
2025-07-08LSTM 輸出不確定的成因、影響與應對策略? 長短期記憶網絡(LSTM)作為循環(huán)神經網絡(RNN)的一種變體,憑借獨特的門控機制,在 ...
2025-07-07統計學方法在市場調研數據中的深度應用? 市場調研是企業(yè)洞察市場動態(tài)、了解消費者需求的重要途徑,而統計學方法則是市場調研數 ...
2025-07-07CDA數據分析師證書考試全攻略? 在數字化浪潮席卷全球的當下,數據已成為企業(yè)決策、行業(yè)發(fā)展的核心驅動力,數據分析師也因此成為 ...
2025-07-07剖析 CDA 數據分析師考試題型:解鎖高效備考與答題策略? CDA(Certified Data Analyst)數據分析師考試作為衡量數據專業(yè)能力的 ...
2025-07-04SQL Server 字符串截取轉日期:解鎖數據處理的關鍵技能? 在數據處理與分析工作中,數據格式的規(guī)范性是保證后續(xù)分析準確性的基礎 ...
2025-07-04CDA 數據分析師視角:從數據迷霧中探尋商業(yè)真相? 在數字化浪潮席卷全球的今天,數據已成為企業(yè)決策的核心驅動力,CDA(Certifie ...
2025-07-04CDA 數據分析師:開啟數據職業(yè)發(fā)展新征程? ? 在數據成為核心生產要素的今天,數據分析師的職業(yè)價值愈發(fā)凸顯。CDA(Certified D ...
2025-07-03