인자분석

인자분석 또는 요인분석

인자 또는 요인이라고 불리는 잠재적으로 적은 숫자의 관찰되지 않은 변수들로, 관찰된 서로 상관인 변수들 사이에서의 분산을 설명하기 위한 통계학적 방법

공통인자 (Common Factor) : $f_{ik}$
인자적재 (Factor Loading) : $\delta_{jk}$
공통성(Communality) : $h^2_j=\sum_{k=1}^{m} \delta_{jk}^2$     ->   공통분산(common variance)
특정분산 (Specifi Variance) : $\Phi_j$ , uniquenesses

R코드

##### 인자분석 #####
D=diag(c(1,2,4,2))
rr=matrix(c(1,0.6,0.1,-0.1,
            0.6,1,0.2,-0.2,
            0.1,0.2,1,0.7,
            -0.1,-0.2,0.7,1),nrow=4, byrow=T)
Sigma=D^(1/2)%*%rr%*%D^(1/2)

library(MASS)
nn=10000
mu=rep(0,dim(rr)[1])
set.seed(2023)
d1=scale(mvrnorm(nn,mu,Sigma))
R=cov(d1)

eeR=eigen(R)

U=eeR$vectors
L=diag(eeR$values)

R
U%*%L%*%t(U)

m=2
Um=U[,1:m]
Lm=L[1:m,1:m]
Rm=Um%*%Lm%*%t(Lm)

round(R-Rm,3)

Umr=U[,-(1:m)]
Lmr=L[-(1:m),-(1:m)]
Rmr=Umr%*%Lmr%*%t(Lmr)

eeR2=eigen(R-diag(diag(Rmr)))
eeR2
11=eeR2$values
cumsum(11)/sum(11)

psi=diag(diag(Rmr))
D=U%*%sqrt(L)

dim(t(D))
dim(solve(psi))

round(solve(t(D)%*%solve(psi)%*%D),3)
f1=solve(t(D)%*%solve(psi)%*%D)%*%t(D)%*%solve(psi)%*%d1[1,]
f1
f=d1%*%solve(psi)%*%D%*%solve(t(D)%*%solve(psi)%*%D)
f[1,]

plot(f[,1],f[,2],pch=19,col='grey')

###
pp=10
set.seed(2022)
Sigma=matrix(rWishart(1,pp,diag(pp)),pp)
R=diag(1/sqrt(diag(Sigma)))%*%Sigma%*%diag(1/sqrt(diag(Sigma)))

nn=10000
mu=rep(0,pp)
set.seed(2022)
d1=mvrnorm(nn,mu,Sigma)
d1=data.frame(d1)
names(d1)=paste('X',1:pp,sep="")
Sigma
cov(d1)

eeS=eigen(cor(d1))
round(cumsum(eeS$values/sum(eeS$values)*100),2)

d1.pca <- princomp(d1,cor=T)
summary(d1.pca)

mm=4
d1.fa=factanal(d1,mm)
loadings(d1.fa)

round(cor(d1),2)
D1=d1.fa$loadings[1:pp,1:mm]
P1=diag(d1.fa$uniquenesses)
round(D1%*%t(D1)+P1,2)

##KMO표본적합성 측도
###overall MSA클수록 요인분석 필요 (0.5기준)
head(d1)
dim(d1)
cov(d1)

#install.packages('psych')
library(psych)
R=cor(d1)
KMO(d1)

s1<-diag(sqrt(1/diag(solve(R))))
q1=s1%*%solve(R)%*%s1
sur2=sum((R-diag(diag(R)))^2)
suq2=sum((q1-diag(diag(q1)))^2)
sur2/(sur2+suq2)

##Bartlett의 구형성 검정
#H_0이 기각인 경유 요인분석 진행
library(psych)
cortest.bartlett(R,nn)

tmp=mvrorm(nn,mu,diag(1,length(mu)))
cortest.bartlett(cor(tmp),nn)

##Wishart 분포
##Hotelling T^2
library(MASS)
pp=8
set.seed(2022)
Sigma=matrix(rWishart(1,pp,diag(pp)),pp)
R=diag(1/sqrt(diag(Sigma)))%*%Sigma%*%diag(1/sqrt(diag(Sigma))) #표본상관행렬

nn=10000
mu=rep(0,pp)
set.seed(2022)
d1=mvrnorm(nn,mu,Sigma)
d1=data.frame(d1)
names(d1)=paste('X',1:pp,sep="")
Sigma

mm=3 #공통인자 수
d1.fa=factanal(scale(d1),mm)
#표준화 
d1.fa=factanal(scale(d1),4)
d1.fa$loadings

d1.fa=factanal(d1,4)
d1.fa$loadings


d1.fa=factanal(scale(d1),3)
d1.fa$loadings
round(d1.fa$loadings[1:pp,],3)
round(d1.fa$uniquenesses,3)
round(1-apply(d1.fa$loadings^2,1,sum),3)

ss1=apply(d1.fa$loadings^2,1,sum) 
round(ss1,3)
round(ss1/(sum(ss1)+sum(d1.fa$uniquenesses)),3) #원래자료의 각 요인의 비율
  #proportion var와 같음

#residual matrix
round(cor(d1),2) #표본상관행렬
D1=d1.fa$loadings[1:pp,1:mm]#인자적재행렬
P1=diag(d1.fa$uniquenesses)#특정분산행렬
round(D1%*%t(D1)+P1,2) #cov(Y)
round(cor(d1)-(D1%*%t(D1)+P1),2)

round(diag(t(D1)%*%D1)/sum(diag(D1%*%t(D1)+P1)),3)

D1
biplot(x=D1,y=D1)


##Rotation
d1.fa.1=factanal(d1,mm,rotation = 'none')
d1.fa.2=factanal(d1,mm,rotation = 'varimax')
d1.fa.3=factanal(d1,mm,rotation = 'promax')
d1.fa.1
d1.fa.2
d1.fa.3

#install.packages("GPArotation")
library(GPArotation)
d1.fa.4=factanal(d1,mm,rotation='quartimax')
d1.fa.5=factanal(d1,mm,rotation='oblimin')



dd=d1.fa.1
dd$lodings[1:8,]
dd2=matrix(as.character(dd),8)
dd2[abs(dd)<0.5]=" "
as.data.frame(dd2)
#해석이 용이해짐

###
apply(loadings(d1.fa.1)^2,2,var)
apply(loadings(d1.fa.2)^2,2,var)
apply(loadings(d1.fa.3)^2,2,var)
apply(loadings(d1.fa.4)^2,2,var)
apply(loadings(d1.fa.5)^2,2,var)
#분산의 합 통해서 변동성 확인
#공통요인 통해서 큰 변화있고 어떤 요인으로 설면되나 확인하고 싶음
#결과상 rotation있을 때 보다 분산이 커짐

apply(loadings(d1.fa.1)^2,1,var)
apply(loadings(d1.fa.2)^2,1,var)
apply(loadings(d1.fa.3)^2,1,var)
apply(loadings(d1.fa.4)^2,1,var)
apply(loadings(d1.fa.5)^2,1,var)
#어떤 변수가 요인들이 몰려있는지 확인
#수 작으면 골고루 요인들이 들어감


par(mfrow=c(2,3))
D1=d1.fa.1$loadings[1:pp,]
biplot(x=D1,y=D1,cex=2)
D1=d1.fa.2$loadings[1:pp,]
biplot(x=D1,y=D1)
D1=d1.fa.3$loadings[1:pp,]
biplot(x=D1,y=D1)
D1=d1.fa.4$loadings[1:pp,]
biplot(x=D1,y=D1)
D1=d1.fa.5$loadings[1:pp,]
biplot(x=D1,y=D1)

par(mfrow=c(2,3))
D1=d1.fa.1$loadings[1:pp,2:3]
biplot(x=D1,y=D1,cex=2)
D1=d1.fa.2$loadings[1:pp,2:3]
biplot(x=D1,y=D1)
D1=d1.fa.3$loadings[1:pp,2:3]
biplot(x=D1,y=D1)
D1=d1.fa.4$loadings[1:pp,2:3]
biplot(x=D1,y=D1)
D1=d1.fa.5$loadings[1:pp,2:3]
biplot(x=D1,y=D1)


#길이는 가중치 제곱의 합
#짧으면 설명을 잘 안함 길수록 설명이 많아지는 것것
###

ro=varimax(loadings(d1.fa.1));ro
round(ro$rotmat%*%t(ro$rotmat),3)
loadings(d1.fa.1)%*%ro$rotmat

d1.fa.2
apply(loadings(d1.fa.2)^2,2,sum)/sum(apply(loadings(d1.fa.2)^2,2,sum)+sum(d1.fa.2$uniquenesses))
apply(loadings(d1.fa.2)^2,2,sum)+d1.fa.2$uniquenesses

d1.fa.3
apply(loadings(d1.fa.3)^2,2,sum)/sum(apply(loadings(d1.fa.3)^2,2,sum)+sum(d1.fa.3$uniquenesses))
apply(loadings(d1.fa.3)^2,2,sum)+d1.fa.3$uniquenesses

###
library(MASS)
library(Matrix)
pp=10
set.seed(2022)
AR=function(rhorho,pppp){
  rrr=(1:pppp)
  ttt=matrix(rrr,nrow=pppp,ncol=pppp,byrow=T)
  return(rhorho^abs(ttt-rrr))
}

Sigma=bdiag(matrix(rWishart(1,pp,AR(0.9,pp)),pp))

###
nn=1000
mu=rep(0,pp)
set.seed(2022)
d1=mvrnorm(nn,mu,Sigma)
d1=data.frame(scale(d1))
name(d1)=paste('X',1:pp,sep='')

mm=3
d1.fa=factanal(d1[,-4],mm,rotation = 'varimax')
print(d1.fa,digits=2,cutoff= .65,sort=TRUE)

인자점수

1. 가중최소제곱법

   $Y_i=\Delta f_i+e_i\\ cov(f_i)=I, cov(e_i)=\Phi\\ f_i=(\hat{\Delta'} \Phi^{-1}\hat{\Delta})^{-1}\hat{\Delta'} \Phi^{-1}Y_i$

2. 회귀법
   정규성이 만족될 때 사용할 수 있음

   $f$

3. Rotation 반영한 인자점수

   $Y_i=(\Delta T)(T'f_i)+ e_i\\ =\Delta_* f_{ik} + e_i\\ \hat{f_{ik}}=(\hat{\Delta'_*} \hat{\Phi^{-1}}\hat{\Delta_*})^{-1}\hat{\Delta'_*} \hat{\Phi^{-1}}Y_i\\ =(T'\hat{\Delta'} \hat{\Phi^{-1}}\hat{\Delta}T)^{-1}T\hat{\Delta'} \hat{\Phi^{-1}}Y_i$

   만약 정규성이 만족된다면

   $\hat{f_{ik}}=\Delta_*'R^{-1}Y_i=T'\Delta'R^{-1}Y_i$