인자분석 : 코드 (1)

요인분석

자료행렬 Y가 요인벡터 f에 의해 설명된다고 생각.
그러나 이 식에서 f를 직접 구할 수 없으므로, f와 Y의 연관성을 나타내는 $\Delta$ 를 먼저 구하기.
따라서 Y의 공분산을 분해하고, 설명 불가능한 부분을 빼 $\Delta$ 를 구해낸다.
이를 이용해 요인벡터 f를 구하며, 이때 나오는 숫자는 각 요인 간의 분석을 진행

ex. 요인 간의 분류, 요인들을 통한 예측 등등…

요인분석 - 패키지 없이

자료 만들기 (1)

D = diag(c(1, 2, 0.4, 2))
data = 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)
rr = matrix(data, nrow = 4, byrow = T)
Sigma = D^{1/2} %*% rr %*% D^{1/2}

n = 10000
mu = rep(0, dim(rr)[1])
set.seed(2023)
d1 = scale(mvrnorm(n, mu, Sigma))

\begin{aligned} Y &= \Delta \cdot f_i + e_i \\ cov(Y) &= \Delta\Delta^{\prime} + \Psi \end{aligned}

$Y$ : 표준화된 자료 = 코드에서의 d1
$\Psi$ : 특정분산행렬 (uniqueness)

# cov(Y) 구해보기
R = cov(d1)

# cov(Y)를 고유값 분해
eeR = eigen(R)
U = eeR$vectors
L = diag(eeR$values)

자료 쪼개기

앞서서 우리는 자료의 공분산인 R을 만들었다.

그런데 요인분석의 식은 $\Delta\Delta^{\prime} + \Psi$ 이므로, R을 두 개로 쪼개어 보자.

R = U\Lambda U^{\prime} = \begin{bmatrix} U_1 & \cdots & U_4 \end{bmatrix} \begin{bmatrix} \delta_1 & \\ & \delta_2 & & \\ & & \delta_3 & \\ & & & \delta_4 \\ \end{bmatrix} \begin{bmatrix} U_1 \\ \vdots \\ U_4 \end{bmatrix}

= \begin{bmatrix} U_1 & U_2 \end{bmatrix} \begin{bmatrix} \delta_1 & \\ & \delta_2 \\ \end{bmatrix} \begin{bmatrix} U_1 \\ U_2 \end{bmatrix} + \begin{bmatrix} U_3 & U_4 \end{bmatrix} \begin{bmatrix} \delta_3 & \\ & \delta_4 \\ \end{bmatrix} \begin{bmatrix} U_3 \\ U_4 \end{bmatrix}

# R의 1st, 2nd 자료
m = 2
Um = U[, 1:m]
Lm = L[1:m, 1:m]
Rm = Um %*% Lm %*% t(Um); Rm
round(R - Rm, 3)

## R의 3rd, 4th 자료
Umr = U[, -(1:m)]
Lmr = L[-(1:m), -(1:m)]
Rmr = Umr %*% Lmr %*% t(Umr)
round(Rmr, 3)
diag(diag(Rmr))

R = \begin{bmatrix} U_1 & U_2 \end{bmatrix} \begin{bmatrix} \delta_1 & \\ & \delta_2 \\ \end{bmatrix} \begin{bmatrix} U_1 \\ U_2 \end{bmatrix} + \begin{bmatrix} U_3 & U_4 \end{bmatrix} \begin{bmatrix} \delta_3 & \\ & \delta_4 \\ \end{bmatrix} \begin{bmatrix} U_3 \\ U_4 \end{bmatrix} = \Delta\Delta^{\prime} + \Psi

원 자료의 전체 공분산 = cov(d1) = R
덧셈식의 앞파트(1st, 2nd 자료) = 주요요인들에 의해 설명된 공분산 = Rm
덧셈식의 뒷파트(3rd, 4th 자료) = 주요요인들에 의해 설명되지 못한 공분산 = Rmr
- 우리는 이것이 굉장히 적기를 바란다. 따라서 덧셈식의 앞파트만으로 전체 R이 설명가능하기를 바란다.

식 변형 후

만약 덧셈식의 뒷파트가 대각행렬이 아니라면, 그것의 주대각원소만 남겨둔다.
- 왜? 덧셈식의 뒷파트 = 특정분산행렬 = 항상 주대각원소만 존재.
- 주대각원소가 아닌 다른 부분이 있다는 것은
- diag(diag(A)) = A의 주대각원소를 뺀 나머지가 전부 0인 행렬
전체 공분산행렬 R에서 이 부분을 뺀 것을 $\Delta\Delta^{\prime}$ 이라고 생각할 수 있을 것이다.

R - diag\Bigg( diag \bigg( \begin{bmatrix} U_3 & U_4 \end{bmatrix} \begin{bmatrix} \delta_3 & \\ & \delta_4 \\ \end{bmatrix} \begin{bmatrix} U_3 \\ U_4 \end{bmatrix} \bigg) \Bigg) \approx \Delta\Delta^{\prime}

~=~ U_\star \Lambda_\star {U_\star}^{\prime} ~=~ U_\star~ {\Lambda_\star}^{1/2} \cdot {\Lambda_\star}^{1/2}~ {U_\star}^{\prime} ~=~ \Delta\Delta^{\prime}

이를 코드로 쓰면 다음과 같을 것이다. 그 결과를 dd.T라고 하자.

dd.T의 변수 중 일부를 다시 골라내보자. (=변수축소)

dd.T = R - diag(diag(Rmr))

## 고유값 분해 후 다시 변수축소
eeR_2 = eigen(ddT) 
## 설명비율로 변수 고르기
l1 = eeR_2$values
cumsum(l1) / sum(l1)

## 2개를 고르는 것으로 보고 근사값 도출
U = eeR_2$vectors[, 1:2]
L = diag(eeR_2$values[1:2])
approx.dd.T = U %*% L %*% t(U)

“변수축소한 값(= approx.dd.T) + Rmr의 대각행렬” 의 합을 구해보자.

이 합이 원래의 R과 얼마나 비슷할지를 한번 생각해보자. 어느 정도 설명을 하고 있음을 알 수 있다.

# R과 비교해보기
round(approx.dd.T + diag(diag(Rmr)), 3)
round(R, 3)

> round(approx.dd.T + diag(diag(Rmr)), 3)
       [,1]   [,2]  [,3]   [,4]
[1,]  0.871  0.690 0.165 -0.136
[2,]  0.690  0.926 0.163 -0.163
[3,]  0.165  0.163 0.944  0.764
[4,] -0.136 -0.163 0.764  0.947
> round(R, 3)
       [,1]   [,2]  [,3]   [,4]
[1,]  1.000  0.596 0.106 -0.079
[2,]  0.596  1.000 0.194 -0.191
[3,]  0.106  0.194 1.000  0.709
[4,] -0.079 -0.191 0.709  1.000

델타와 psi 도출

그래서, 이를 기반으로 특정분산행렬 $\Psi$ = psi와 $\Delta$ = D 를 만들었다.

U_\star~ {\Lambda_\star}^{1/2} \cdot {\Lambda_\star}^{1/2}~ {U_\star}^{\prime} ~=~ \Delta\Delta^{\prime}

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

요인점수 계산

오차를 반영할 수 있는 WLS(가중 최소제곱법)을 이용한다.

\underline\beta = {(X^{\prime}~\Sigma^{-1}~X)}^{-1} ~\cdot X^{\prime}\cdot\Sigma^{-1}\cdot Y

f ={(\Delta^{\prime}~\Psi^{-1}~\Delta)}^{-1} \cdot \Delta^{\prime} \cdot \Psi^{-1}\cdot Y

이것을 통해 ${(\Delta^{\prime}~\Psi^{-1}~\Delta)}^{-1}$ 를 도출해보면 다음과 같다.

> round(solve(t(D) %*% solve(psi) %*% D), 3)
      [,1]  [,2]
[1,] 0.069 0.005
[2,] 0.005 0.133

여기서는 첫번째 인자점수를 뽑아보도록 하자.
역행렬 계산 = solve(), 행렬곱 = %*% .

> f = solve(t(D) %*% solve(psi) %*% D) %*% t(D) %*% solve(psi) %*% d1[1,]
> f
           [,1]
[1,] -1.4442125
[2,] -0.3464697

요인분석 - 패키지로

위샤트분포를 이용한 공분산행렬

위샤트분포 : 표본공분산 행렬의 분포
R에서 위샤트분포로 자료 뽑기 : rWishart()

library(MASS)
pp = 10 # 10명
set.seed(2022)
Sigma = matrix(rWishart(1, pp, diag(pp)), pp)

이것을 기반으로 자료를 생성.

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

만들어낸 자료를 확인해보자.

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

주성분분석 vs 인자분석

위에서 만들어낸 자료를 가지고 주성분분석과 인자분석을 동시에 수행.

> d1.pca = princomp(d1, cor=T)
> summary(d1.pca)
Importance of components:
                          Comp.1    Comp.2    Comp.3
Standard deviation     1.8375260 1.4354835 1.2461378
Proportion of Variance 0.3376502 0.2060613 0.1552859
Cumulative Proportion  0.3376502 0.5437115 0.6989974
                          Comp.4     Comp.5     Comp.6
Standard deviation     1.0933010 0.91508227 0.78134893
Proportion of Variance 0.1195307 0.08373756 0.06105062
Cumulative Proportion  0.8185281 0.90226570 0.96331631
                           Comp.7     Comp.8     Comp.9
Standard deviation     0.44793512 0.34993761 0.18779190
Proportion of Variance 0.02006459 0.01224563 0.00352658
Cumulative Proportion  0.98338090 0.99562653 0.99915311
                          Comp.10
Standard deviation     0.09202663
Proportion of Variance 0.00084689
Cumulative Proportion  1.00000000

주성분 4개의 누적비율이 약 0.81으로, 4개가 적합해 보인다!

이것을 scree plot으로도 확인해보자.

## scree plot 
plot(1:10, d1.pca$sdev ^ 2)

elbow를 찾기 조금 애매해 보인다. (실제로 애매한 자료가 많다.)
꺽은 회귀직선, 70% 설명비율, karlson’s method 등 여려가지 방법을 생각해 볼 수 있을 것 같다.

인자분석

> loadings(d1.fa)

Loadings:
      Factor1 Factor2 Factor3 Factor4
 [1,]  0.203  -0.710   0.448   0.322 
 [2,]  0.944                  -0.306 
 [3,] -0.455   0.827  -0.269  -0.151 
 [4,]          0.169   0.970  -0.158 
 [5,]  0.805           0.318         
 [6,]  0.578           0.379   0.548 
 [7,] -0.345   0.208   0.105         
 [8,] -0.222  -0.106  -0.110   0.894 
 [9,]          0.102  -0.355         
[10,]  0.115   0.906   0.269   0.116 


               Factor1 Factor2 Factor3 Factor4
SS loadings      2.312   2.114   1.687   1.363
Proportion Var   0.231   0.211   0.169   0.136
Cumulative Var   0.231   0.443   0.611   0.748

Loadings → 인자적재행렬 \Delta에 해당
- 4개의 요인이 주 영향력을 미치는 인자가 무엇인지
- blank → 임의로 조정이 가능하다.

인자분석의 필요성

만약 $cov(Y) = R$ 이 대각행렬이면, 인자적재행렬 $\Delta$ 가 필요 없음.

왜? 자료간의 상관관계가 없다는 뜻이기 때문에!

cov(Y) = \begin{bmatrix} a & 0 & \cdots & 0 \\ 0 & b & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & z \\ \end{bmatrix} = \Delta\Delta^{\prime} + \Psi

따라서 f로 설명이 되지 않기 때문에, 인자분석의 필요가 없어짐.

R이 대각행렬임을 검증해보는 두 가지 방법이 존재함.

kmo 지수

> # install.packages("psych")
> library(psych)
> KMO(d1)
Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = d1)
Overall MSA = 0.3
MSA for each item = 
 [1] 0.36 0.26 0.33 0.23 0.35 0.52 0.30 0.17 0.16 0.24

Overall Msa가 0.5 이상이면 요인분석이 필요한 것으로 해석하면 됨.

이것의 원리는 다음과 같다.

s1 = diag(sqrt(1 / diag(solve(R))))
q1 = s1 %*% solve(R) %*% s1

sum_r2 = sum((R - diag(diag(R))) ^ 2)
sum_q2 = sum((q1 - diag(diag(q1))) ^ 2)
kmo_value = sum_r2 / (sum_r2 + sum_q2)
> kmo_value
[1] 0.3005988

s1 : 상관행렬 R의 역행렬을 계산하고, 거기서 대각성분을 추출. 그 대각성분들의 역수에 제곱근을 씌운다.
q1 : (s1) x (r의 역행렬) x (s1). 원래의 상관행렬에 변수 간의 부분 상관행렬을 조정
sum_r2 : R에서 대각선을 제외한 모든 부분의 제곱을 씌우고, 그것들을 전부 더한 것
sum_q2 : q1에서 대각선을 제외한 모든 부분의 제곱을 씌옥, 그것들을 전부 더한것
- 부분 상관 계수의 제곱의 합으로, 각 변수가 다른 변수의 영향을 배제했을 때 남는 상관관계의 강도
kmo_value : 전체 상관관계의 제곱합(sum_r2) 대비 전체 상관관계의 제곱합와 부분 상관관계의 제곱합 (sum_r2 + sum_q2)을 구한 값.

Bartelett의 구형성 검정

요인분석의 필요성을 측정하는 가설검정. 가능도비 검정으로 도출한 카이제곱검정의 일종
$H_0 : R = I$ 라고 놓고, 이것이 기각되면 요인분석이 필요하다고 생각.

> cortest.bartlett(R, mm)
$chisq
[1] -11.6333

$p.value
[1] 1

$df
[1] 45

p-value를 보고 이것이 유의수준보다 높다면, 요인분석을 하기에 부적합하다고 해석하면 됨.
귀무가설 하에서의 모수의 갯수 vs 대립가설 하에서의 모수의 갯수를 이용해서 가능도비를 추출