高校生向けの主成分分析教材.
参考テキスト 統計科学研究所 様.　Rによる主成分分析.

1 分散共分散行列.

1変量確率変数の分散はイメージしやすい。分散が大きければ広がりが大きい。
２変量確率変数の分散共分散行列から分布の状態をイメージできるだろうか。
固有値・固有ベクトルとの関係も含めてシミュレーションしてみたい。
2変量正規分布に対して主成分分析を適用する。

共分散とはX,Yの一致度である。
共分散行列とは？(定義、意味、分布の推定).

1.1 相関のある2変量

簡単のため、2変量正規乱数を使用する。
算出するのは、

  • 分散共分散行列.
  • 固有値.
  • 固有ベクトル.

散布図と固有ベクトルの関係をつかもう。

library(MASS)
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.7, 0.7, 1), 2, 2)
data <- mvrnorm(50, mu, Sigma)
d <- data.frame(X=data[,1],Y=data[,2])
cov(d)#分散共分散行列
##           X         Y
## X 1.1747491 0.6327685
## Y 0.6327685 0.9769903
eigen_d <- eigen(cov(d))#固有値
v <- eigen_d$vectors
e <- eigen_d$values
plot(d,asp=1)
arrows(x0=0, y0=0, x1=e[1]*v[1,1],  y1=e[1]*v[2,1], code=2, col="red")
arrows(x0=0, y0=0, x1=v[1,2],  y1=v[2,2], code=2, col="blue")

分布は、
　固有ベクトル (1,1) 方向に 1+r 倍.
　固有ベクトル (1,-1) 方向に 1-r 倍.
相関係数rが大きくなれば楕円の短軸は圧縮される。

1.2 主成分分析.

library("psych")
#(describe(d))#基本統計情報
par(family = "HiraKakuProN-W3") #日本語フォントの指定
result <- prcomp(d, scale=T) #scale=Tは相関行列から分析
summary(result)
## Importance of components:
##                           PC1    PC2
## Standard deviation     1.2612 0.6398
## Proportion of Variance 0.7953 0.2047
## Cumulative Proportion  0.7953 1.0000
biplot(result)

biplot は主成分得点と主成分負荷量を同時にプロットしたもので、座標軸に２種類の目盛がある。便利であるがテキストをグラフに記入するため表示にズレがあるので、次のように plotで描くこともできる。

主成分負荷量は主成分との相関係数であり、rotationの行をベクトル表示したものである。主成分が元の変数に影響する度合いを図から読み取ることができる。
因子負荷量をプロットすることにより、主成分に寄与している因子を視覚的に捉えることができる。

plot(result$x[,1],result$x[,2],type="n",asp=1) 
text(result$x[,1],result$x[,2],1:50)

主成分の計算を再現する。
直交行列result$rotation の１列目が第１主成分を求める式の係数。
dを標準化したbの１行目が標準化された元のデータ。
内積を計算すると第１主成分となる。

result$x[1,1:2]#主成分
##        PC1        PC2 
##  0.9020726 -0.9641003
a<-result$rotation;a
##         PC1        PC2
## X 0.7071068  0.7071068
## Y 0.7071068 -0.7071068
b<-scale(d)
a[1:2,1]
##         X         Y 
## 0.7071068 0.7071068
b[1,1:2]
##           X           Y 
## -0.04386021  1.31958347
a[1:2,1]*b[1,1:2]
##           X           Y 
## -0.03101385  0.93308642
sum(a[1:2,1]*b[1,1:2])
## [1] 0.9020726

1.3 因子負荷量.

MARGIN=2 は行の操作を表す。（MARGIN=1は列操作）

result$sdev
## [1] 1.2612081 0.6398079
fc.l <- sweep(result$rotation, MARGIN=2, result$sdev, FUN="*") 
result$rotation
##         PC1        PC2
## X 0.7071068  0.7071068
## Y 0.7071068 -0.7071068
z<-fc.l#主成分負荷行列
z
##         PC1        PC2
## X 0.8918088  0.4524125
## Y 0.8918088 -0.4524125
z[,1]#第１主成分の主成分負荷量
##         X         Y 
## 0.8918088 0.8918088
sum(z[,1]^2)#大きさは第１固有値
## [1] 1.590646
par(family = "HiraKakuProN-W3") #日本語フォントの指定
subject <- c("X", "Y") ;subject
## [1] "X" "Y"
plot(fc.l[,1], pch=subject, ylim=c(-1,1), main="PC1")

plot(fc.l[,2], pch=subject, ylim=c(-1,1), main="PC2")

#主成分との相関係数

1.4 バリマックス回転の主成分分析.

library(psych)
dat <- d
fit <- principal(dat, nfactors=2, rotate=TRUE)#バリマックス回転
## Specified rotation not found, rotate='none' used
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 2, rotate = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1   PC2 h2 u2 com
## X 0.89 -0.45  1  0 1.5
## Y 0.89  0.45  1  0 1.5
## 
##                        PC1  PC2
## SS loadings           1.59 0.41
## Proportion Var        0.80 0.20
## Cumulative Var        0.80 1.00
## Proportion Explained  0.80 0.20
## Cumulative Proportion 0.80 1.00
## 
## Mean item complexity =  1.5
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0 
##  with the empirical chi square  0  with prob <  NA 
## 
## Fit based upon off diagonal values = 1
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 2, rotate = TRUE)
## 
## Test of the hypothesis that 2 factors are sufficient.
## The degrees of freedom for the model is -2  and the objective function was  0 
## The number of observations was  50  with Chi Square =  0  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0
biplot(fit)

library(psych)
dat <- d
fit <- principal(dat, nfactors=1, rotate=TRUE)#バリマックス回転
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 1, rotate = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1  h2  u2 com
## X 0.89 0.8 0.2   1
## Y 0.89 0.8 0.2   1
## 
##                 PC1
## SS loadings    1.59
## Proportion Var 0.80
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 component is sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.2 
##  with the empirical chi square  4.19  with prob <  NA 
## 
## Fit based upon off diagonal values = 0.88
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 1, rotate = TRUE)
## 
## Test of the hypothesis that 1 factor is sufficient.
## The degrees of freedom for the model is -1  and the objective function was  0.31 
## The number of observations was  50  with Chi Square =  14.7  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0.2
biplot(fit)

2 3変量.

#3D表示
options(rgl.printRglwidget = TRUE) 
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)

library(MASS)
mu <- c(0,0,0)
Sigma <- matrix(c(1, 0.7,0.6,0.7,1,0.5,0.6,0.5, 1), 3, 3)
data <- mvrnorm(100, mu, Sigma)
d <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
cov(d)#分散共分散行列
##           X         Y         Z
## X 0.8399898 0.5724266 0.5052255
## Y 0.5724266 0.7513178 0.3403845
## Z 0.5052255 0.3403845 0.8662971
attach(d)
#plot3d(X, Y, Z, type = "s", col = "blue", size = 1)

eigen_d <- eigen(cov(d))#固有値
v <- eigen_d$vectors;v[1,1]
## [1] -0.632473
e <- eigen_d$values;e[1]
## [1] 1.772291
plot3d(X, Y, Z, type = "s", col = "gray", size = 1)
arrow3d(p0 = c(0, 0, 0), p1 = c(e[1]*v[1,1],e[1]*v[2,1],e[1]*v[3,1]), col = "red")
arrow3d(p0 = c(0, 0, 0), p1 = c(2*e[2]*v[1,2],2*e[2]*v[2,2],2*e[2]*v[3,2]), col = "blue")
arrow3d(p0 = c(0, 0, 0), p1 = c(3*e[3]*v[1,3],3*e[3]*v[2,3],3*e[3]*v[3,3]), col = "green")
library("psych")
#(describe(d))#基本統計情報
par(family = "HiraKakuProN-W3") #日本語フォントの指定
result <- prcomp(d, scale=T) #scale=Tは相関行列から分析
summary(result)
## Importance of components:
##                           PC1    PC2     PC3
## Standard deviation     1.4714 0.7689 0.49383
## Proportion of Variance 0.7217 0.1971 0.08129
## Cumulative Proportion  0.7217 0.9187 1.00000
biplot(result)

plot(result$x[,1],result$x[,2],type="n",asp=1) 
text(result$x[,1],result$x[,2],1:100) 

library(psych)
dat <- d
fit <- principal(dat, nfactors=2, rotate=TRUE)#バリマックス回転
## Specified rotation not found, rotate='none' used
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 2, rotate = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1   PC2   h2    u2 com
## X 0.92 -0.11 0.86 0.144 1.0
## Y 0.85 -0.45 0.92 0.081 1.5
## Z 0.77  0.62 0.98 0.019 1.9
## 
##                        PC1  PC2
## SS loadings           2.16 0.59
## Proportion Var        0.72 0.20
## Cumulative Var        0.72 0.92
## Proportion Explained  0.79 0.21
## Cumulative Proportion 0.79 1.00
## 
## Mean item complexity =  1.5
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.07 
##  with the empirical chi square  3.19  with prob <  NA 
## 
## Fit based upon off diagonal values = 0.98
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 2, rotate = TRUE)
## 
## Test of the hypothesis that 2 factors are sufficient.
## The degrees of freedom for the model is -2  and the objective function was  0.41 
## The number of observations was  100  with Chi Square =  38.87  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0.07
biplot(fit)

library(psych)
dat <- d
fit <- principal(dat, nfactors=3, rotate=TRUE)#バリマックス回転
## Specified rotation not found, rotate='none' used
fit
## Principal Components Analysis
## Call: principal(r = dat, nfactors = 3, rotate = TRUE)
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1   PC2   PC3 h2       u2 com
## X 0.92 -0.11 -0.38  1  2.2e-16 1.4
## Y 0.85 -0.45  0.28  1 -1.2e-15 1.8
## Z 0.77  0.62  0.14  1 -2.2e-16 2.0
## 
##                        PC1  PC2  PC3
## SS loadings           2.16 0.59 0.24
## Proportion Var        0.72 0.20 0.08
## Cumulative Var        0.72 0.92 1.00
## Proportion Explained  0.72 0.20 0.08
## Cumulative Proportion 0.72 0.92 1.00
## 
## Mean item complexity =  1.7
## Test of the hypothesis that 3 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0 
##  with the empirical chi square  0  with prob <  NA 
## 
## Fit based upon off diagonal values = 1
summary(fit)
## 
## Factor analysis with Call: principal(r = dat, nfactors = 3, rotate = TRUE)
## 
## Test of the hypothesis that 3 factors are sufficient.
## The degrees of freedom for the model is -3  and the objective function was  0 
## The number of observations was  100  with Chi Square =  0  with prob <  NA 
## 
## The root mean square of the residuals (RMSA) is  0
biplot(fit)