위샤트분포
자료 만들기
r
## 위샤트분포로 공분산행렬 제작
library(MASS)
pp = 8
set.seed(2022)
Sigma = matrix(rWishart(1, pp, diag(pp)), pp)
r
## 표본상관행렬
R = diag(1 / sqrt(diag(Sigma))) %*% Sigma %*% diag(1/sqrt(diag(Sigma)))
rWishart()
1: 표본의 갯수pp: 차원을 의미 (각 subject별 8개의 자료)diag(pp): 초깃값을 의미. (8 x 8의 단위행렬)
R : 표본상관행렬
- Sigma를 표준편차로 2번 나눠줬다고 생각하면 됨.
r
## 공분산행렬 Sigma를 확인
round(Sigma, 2)
## 양정치행렬 확인
if (det(Sigma) > 0) paste("양정치행렬") else paste("양정치행렬이 아님")
## 상관행렬을 확인
> round(R, 2)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1.00 0.21 -0.89 0.29 0.41 0.46 -0.15 -0.12
[2,] 0.21 1.00 -0.51 0.15 0.84 0.48 -0.40 -0.48
[3,] -0.89 -0.51 1.00 -0.17 -0.55 -0.49 0.35 0.45
[4,] 0.29 0.15 -0.17 1.00 0.36 0.39 0.16 0.69
[5,] 0.41 0.84 -0.55 0.36 1.00 0.76 -0.38 -0.19
[6,] 0.46 0.48 -0.49 0.39 0.76 1.00 -0.32 0.06
[7,] -0.15 -0.40 0.35 0.16 -0.38 -0.32 1.00 0.35
[8,] -0.12 -0.48 0.45 0.69 -0.19 0.06 0.35 1.00
> 1st, 3rd 간에는 강한 음의 상관관계가 존재.
r
## 자료 만들기
n = 10000
mu = rep(0, pp)
set.seed(2022)
d1 = mvrnorm(n, mu, Sigma)
d1 = data.frame(d1)
names(d1) = paste('X', 1:pp, sep='')
요인분석 예시 (1)
공통요인 3개
r
> mm = 3
> d1.fa = factanal(scale(d1), mm)
> d1.fa
Call:
factanal(x = scale(d1), factors = mm)
Uniquenesses:
X1 X2 X3 X4 X5 X6 X7 X8
0.066 0.005 0.005 0.091 0.191 0.519 0.790 0.005
Loadings:
Factor1 Factor2 Factor3
X1 0.126 0.952 0.109
X2 0.982 0.105 -0.141
X3 -0.394 -0.903 0.157
X4 0.269 0.234 0.884
X5 0.839 0.283 0.156
X6 0.488 0.388 0.304
X7 -0.352 -0.180 0.230
X8 -0.325 -0.181 0.926
- mm = 공통인자의 갯수.
- Uniquenesses : \Psi에 해당하는 부분. 1에 가까워질수록 요인 f가 설명하는 부분이 적어진다는 의미.
- 공통인자로 설명되지 못하는 부분의 분산.
- 7번째에 해당하는 부분(
X7)은 현재 공통인자로 잘 설명되지 못한다.
- 그런데, X7에 해당하는 부분이 중요하다면?
- 요인의 갯수를 늘려야 한다.
공통요인 4개
r
> mm = 4
> d1.fa = factanal(scale(d1), mm); d1.fa
Call:
factanal(x = scale(d1), factors = mm)
Uniquenesses:
X1 X2 X3 X4 X5 X6 X7 X8
0.034 0.066 0.005 0.005 0.019 0.289 0.724 0.005
여전히 X7이 잘 설명되지 않고 있다.
공통요인 5개
r
> mm = 5
> d1.fa = factanal(scale(d1), mm)
Error in factanal(scale(d1), mm) : 5 factors are too many for 8 variables
- \Psi가 대각행렬이 되어야 하는데, 그 구조가 무너졌기에 결과가 나오지 않는다.
요인분석 예시 (2)
이번에는 변수를 빼고 진행해보자.
공통요인 3개 + 6, 7번째 변수 빼기
r
mm = 3
d1.fa = factanal(d1[, -c(6,7)], mm); d1.fa
Call:
factanal(x = d1[, -c(6, 7)], factors = mm)
Uniquenesses:
X1 X2 X3 X4 X5 X8
0.066 0.005 0.005 0.089 0.194 0.005
공통요인 3개 + 7번째 변수 빼기
r
## 공통요인 3개 + 7번째 변수 빼기 => x6가 잘 설명이 안된다.
mm = 3
d1.fa = factanal(d1[, -7], mm); d1.fa
Call:
factanal(x = d1[, -7], factors = mm)
Uniquenesses:
X1 X2 X3 X4 X5 X6 X8
0.066 0.005 0.005 0.091 0.192 0.519 0.005
공통요인 4개 + 7번째 변수 빼기
r
## 공통요인 3개 + 7번째 변수 빼기 => 에러 발생
> mm = 4
> d1.fa = factanal(d1[, -7], mm)
Error in factanal(d1[, -7], mm) : 4 factors are too many for 7 variables
r
## 최종결과
mm = 4
d1.fa = factanal(scale(d1), mm); d1.fa
Call:
factanal(x = scale(d1), factors = mm)
Uniquenesses:
X1 X2 X3 X4 X5 X6 X7 X8
0.034 0.066 0.005 0.005 0.019 0.289 0.724 0.005
Loadings:
Factor1 Factor2 Factor3 Factor4
X1 0.942 0.267
X2 0.117 -0.115 0.504 0.808
X3 -0.902 0.123 -0.253 -0.320
X4 0.214 0.939 0.154 0.210
X5 0.227 0.819 0.500
X6 0.316 0.191 0.750 0.112
X7 -0.161 0.301 -0.365 -0.161
X8 -0.214 0.878 -0.422
특정분산 계산하기
특정분산(uniqueness) 산정방법
\begin{aligned}
Y_i &= \Delta f_i + e_i \\
cov(Y_i) &= \Delta \Delta^{\prime} + \Psi \\[10pt]
Y_i &=
\begin{bmatrix}
Y_{i1} \\ \vdots \\ Y_{ip}
\end{bmatrix}_,~
e_i =
\begin{bmatrix}
e_{i1} \\ \vdots \\ e_{ip}
\end{bmatrix}
\\[25pt]
Y_{i1} &= \delta_{11}~f_1 + \delta_{12}~f_2 + \cdots +
\delta_{1m}~f_m+ e_{i1} \\
Var(Y_{i1}) &= {\delta_{11}}^2 + {\delta_{12}}^2 + \cdots + {\delta_{1m}}^2 + \psi_1 = 1 \\
\psi_1 &= 1 - ({\delta_{11}}^2 + {\delta_{12}}^2 + \cdots + {\delta_{1m}}^2)
\end{aligned}
코드로 옮기기
r
# 요인분석
mm = 3
d1.fa = factanal(scale(d1), mm)
## 패키지의 uniqueness
unique1 = round(d1.fa$uniquenesses, 3)
## 직접 구한 uniqueness
d1.fa$loadings
round(d1.fa$loadings[1:pp], 3)
unique2 = 1 - apply(d1.fa$loadings^2, 1, sum
unique2 = round(unique2, 3)
r
# apply(d1.fa$loadings^2, 1, sum) 가 의미하는 것
Loadings:
Factor1 Factor2 Factor3
X1 0.126 0.952 0.109
X2 0.982 0.105 -0.141
X3 -0.394 -0.903 0.157
X4 0.269 0.234 0.884
X5 0.839 0.283 0.156
X6 0.488 0.388 0.304
X7 -0.352 -0.180 0.230
X8 -0.325 -0.181 0.926
unique2 = 1 - ([0.126 0.952 0.109] 각각을 전부 제곱한 뒤 더함)
* 두번째 파리미터의 1 : column 방향을 의미
r
> if (all(unique1 == unique2)) paste("같다") else paste("다르다")
[1] "같다"
분산의 비율 (Proportion Var.)
Proportion Var를 손으로 도출하기.
r
ssl = apply(d1.fa$loadings^2, 2, sum)
-------------------------------------------------
Factor1 Factor2 Factor3
X1 0.126 0.952 0.109
X2 0.982 0.105 -0.141
X3 -0.394 -0.903 0.157
X4 0.269 0.234 0.884
X5 0.839 0.283 0.156
X6 0.488 0.388 0.304
X7 -0.352 -0.180 0.230
X8 -0.325 -0.181 0.926
ssl[1] : [0.126 0.982 ... -0.352 -0.325]을 각각 제곱 후 합친 값.
ssl[2] : [0.952 0.105 ... -0.180 -0.181]을 각각 제곱 ㅎ 합친 값.
...
r
> ssl = round(ssl, 3)
> round(ssl / (sum(ssl) + sum(d1.fa$uniquenesses)), 3)
Factor1 Factor2 Factor3
0.297 0.260 0.233
이것은 아래의 굵은 글씨로 표기한 부분과 같다.
r
> d1.fa$loadings
Loadings: ....
Factor1 Factor2 Factor3
SS loadings 2.379 2.084 1.865
Proportion Var 0.297 0.260 0.233
Cumulative Var 0.297 0.558 0.791
- 공통인자로 전체 자료를 79% 정도 설명 가능하다.
만약 공통인자의 갯수가 4개라면? 좋은 기준이라고 생각은 할 수 있다.
r
> mm = 4
> d1.fa = factanal(scale(d1), mm)
> d1.fa$loadings
Factor1 Factor2 Factor3 Factor4
Cumulative Var 0.248 0.476 0.698 0.857
공통인자로 설명할 수 없는 부분 = 남은 15% = uniqueness
그 부분이 X1, X2, X6, X7 이라고 유추 가능.
r
> d1.fa$uniquenesses
X1 X2 X3 X4 X5
0.03363592 0.06621804 0.00500000 0.00500000 0.01888823
X6 X7 X8
0.28863945 0.72429468 0.00500000
r
## tip : 생활의 지혜
d1.fa$loadings # 행렬 2개 도출
d1.fa$loadings[1:3, ] # 헹렬 1개만 도출. 수치화 가능
잔차행렬
(원래 표본상관행렬) - (공통인자로 만든 표본상관행렬)의 차이를 의미
r
## 원래의 표본상관행렬
mat1 = round(cor(d1), 2)
## 공통인자로 추정한 표본상관행렬
pp = 8; mm = 3
D1 = d1.fa$loadings[1:pp, 1:mm]
P1 = diag(d1.fa$uniquenesses)
mat2 = round(D1 %*% t(D1) + P1, 2)
## 둘의 차이 보기 : 약간의 차이 존재
mat1;mat2
둘의 차를 이용해서 잔차행렬을 구해보자.
r
mat1 = cor(d1)
residual_matrix = round((mat1 - mat2), 2)
residual_matrix
설명이 안되는 변수들과 관련된 부분들이 차이가 크다는 것을 알 수 있음.
인자적재행렬을 이용한 biplot
r
round(D1, 2)
par(mfrow=c(1, 1))
biplot(x = D1, y = D1)
> round(D1, 2)
Factor1 Factor2 Factor3
X1 0.94 0.08 0.27
X2 0.12 -0.11 0.50
X3 -0.90 0.12 -0.25
X4 0.21 0.94 0.15
X5 0.23 0.09 0.82
X6 0.32 0.19 0.75
X7 -0.16 0.30 -0.37
X8 -0.21 0.88 0.00
x1 = (0.94, 0.08), x2 = (0.12, -0.11) 라고 해석할 수 있다.
따라서 화살표의 방향과 길이는 인자적재행렬 \Delta에서 나오는 것이라고 볼 수 있다.
Cov(Y_{i1}, f_{i1}) = cov(\delta_{i1}, f_{i1} + \delta_{im}, f_{im}) = \delta_{i1}
- 화살표의 방향과 평행한 축 = 영향을 많이 미치는 요인
- 길이가 길수록 요인과 관련성이 깊다.