자료 제작
이전 “주성분분석 : 코드(1)”의 4x4 공분산행렬에서 쓴 코드를 사용
r
library(MASS)
D = diag(c(1, 2, 0.4, 0.5))
data = c(1, 0.6, 0.1, 0.3,
0.6, 1, 0.2, 0.4,
0.1, 0.2, 1, 0.5,
0.3, 0.4, 0.5, 1)
rr <- matrix(data, nrow=4, byrow=TRUE)
Sigma = D^{1/2} %*% rr %*% D^{1/2}
mu = rep(0, 4)
n = 100
d1 = mvrnorm(n, mu, Sigma)
d1 = scale(d1)
eeR = eigen(cov(d1))
패키지 (1) - prcomp
prcomp() 안에 표준화된 자료행렬(p1)을 넣으면 객체가 반환되며, 이 안에 여러가지 정보가 있다.
- standard deviation : 고유값의 제곱근.
- Rotation : 고유벡터들을 모은 행렬이라고 생각하면 된다.
- 부호가 반대로 되어도 걱정할 필요가 전혀 없음.
r
> d1.pc = prcomp(d1)
> d1.pc
Standard deviations (1, .., p=3):
[1] 1.5748360 0.7094494 0.6470424
Rotation (n x k) = (3 x 3):
PC1 PC2 PC3
[1,] -0.4312426 0.89865671 0.08028646
[2,] -0.9010111 -0.43358408 0.01356221
[3,] -0.0469987 0.06649039 -0.99667956
> plot(d1.pc)
summary()를 통해 주성분의 분산, 점수, 누적점수를 한번에 볼 수 있다.
r
> summary(d1.pc)
Importance of components:
PC1 PC2 PC3
Standard deviation 1.575 0.7094 0.6470
Proportion of Variance 0.729 0.1479 0.1231
Cumulative Proportion 0.729 0.8769 1.0000
이번에는 d1.pc의 구조를 str()로 보자.
r
> str(d1.pc)
# sdev, rotation, center, scale, x가 있다.
> d1.pc$sdev # 표준편차
[1] 1.5748360 0.7094494 0.6470424
> d1.pc$center # 평균값
[1] 0.03309963 0.02035148 0.05805797
> d1.pc$scale # 표준화
[1] FALSE
표준편차 = 고유값의 제곱근인가?
실제로 그런지 확인해보자.
r
> ## 표준편차 = 고유값의 제곱근인지 확인
> d1.pc$sdev # 표준편차
[1] 1.4213361 1.0536106 0.6983555 0.6180680
> sqrt(eigen(cor(d1))$values) # 고유값의 제곱근
[1] 1.4213361 1.0536106 0.6983555 0.6180680
자료의 축소
설명력이 작은 고유벡터를 제거해보자.
r
# 고유값들을 tmp에 담고, 그것들을 반올림한 것을 tmp2에 담는다.
tmp = eigen(cor(d1))$vectors
tmp2 = round(tmp,3); tmp2
# 고유값이 0.3보다 작으면 (= 설명력이 낮으면) 제거
tmp3 = matrix(as.character(tmp2),4)
tmp3[abs(tmp2 < 0.3)] = ''
as.data.frame(tmp3)
biplot으로 확인
scale에 따라 표준화 여부가 결정, 이때 biplot의 모양이 달라진다.
r
## biplot 그리기
biplot(d1.pc, pch=19, cex=0.6, scale=T, lwd=3)
biplot(d1.pc, pch=19, cex=0.6, scale=F, lwd=3) #
scale = T: 표준화 o -> 원scale = F:표준화 x -> 타원
패키지 (2) - princomp
r
d1.pc2 = princomp(d1)
str(d1.pc2)
loadings 살펴보기
loadings = 고유벡터
r
> loadings(d1.pc2)
Loadings:
Comp.1 Comp.2 Comp.3 Comp.4
[1,] 0.526 0.511 0.680
[2,] 0.573 0.344 -0.255 -0.698
[3,] 0.391 -0.670 -0.596 0.209
[4,] 0.492 -0.415 0.761
Comp.1 Comp.2 Comp.3 Comp.4
SS loadings 1.00 1.00 1.00 1.00
Proportion Var 0.25 0.25 0.25 0.25
Cumulative Var 0.25 0.50 0.75 1.00
객체 보기
아까 prcomp()의 결과와는 살짝 다른 것을 알 수 있음.
r
> d1.pc2
Call:
princomp(x = d1)
Standard deviations:
Comp.1 Comp.2 Comp.3 Comp.4
1.4392711 1.0276119 0.7250234 0.5539435
4 variables and 100 observations.
biplot 그리기
prcomp()로 biplot을 그렸을 때와 달리 방향이 반대로 나올 수 있으나, 그 결과는 같다.
r
biplot(d1.pc2, pch=19, cex=0.6, scale=F, lwd=3)
이상치 찾아내기
biplot을 그리다 보면 뭔가 이상한 부분을 찾을 때가 많다.
이런 경우, 대부분 자료를 이루는 모집단이 두 개 이상인 경우가 많다.
자료 제작 - 모집단 2개 섞기
먼저 공분산행렬과 평균벡터를 만든다.
r
library(MASS)
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=TRUE)
D = diag(c(1, 2, 0.4, 0.2))
Sigma = D^{1/2} %*% rr %*% D^{1/2}
nn = 10000
mm = 20
mu = rep(0, dim(rr)[1])
# mu2 = rep(7, dim(rr)[1])
mu2 = c(7, 0, 1, 2)
rbind()를 사용해 두 집단을 섞는다.
r
n1 = mvrnorm(nn - mm, mu, Sigma)
n2 = mvrnorm(mm, mu2, Sigma)
d1 = rbind(n1, n2)
d1 = scale(d1)
새로운 자료인 d1을 살펴보자.
r
class(d1)
R = cov(d1)
eeR = eigen(R)
cumsum(eeR$values) / sum(eeR$values) * 100
round(eeR$vectors, 2)
이상치 찾아내기
r
P = data.frame(d1 %*% eeR$vectors[, 1:2])
plot(P$X1, P$X2, pch=19, col="grey", cex=0.2,
xlab="1st pc", ylab="2nd pc", bty="l")
