要約

Rの主成分分析は果たして何をしているのか? ちょっと確かめてみました.結論は,以下の通り

  1. データの主成分分析は,各点からの垂線の距離(情報損失量)が最小になるように,固有ベクトルを決める.
  2. 主成分得点の分散の最大化も,これと同じことをしている.
  3. 固有値は,最小化された情報損失量であり,最大化された主成分得点の分散でもある.
  4. Rで主成分分析はprcomp()関数を使う.
  5. 固有値の平方根(主成分得点の標準偏差)$sdevとノルム1で直交する固有ベクトル$rotationが出力される.
  6. 主成分得点(スコア)は$xで出力される.
  7. 通常,主成分分析では,各変数の分散は1に標準化した方がよいが,これは引数scale=Tでできる.
  8. 主成分負荷量は,標準化した変数と,それを使った主成分分析のスコアとの相関係数とするのが一般的.
  9. biplot()関数を使うと,スコアと主成分負荷量という別種の値が1つのグラフにプロットされるが,これらの値は,特異値分解で得られるもので,上記のそれらとは値が異なる.
  10. biplot()関数によるプロット図は,引数scale=に0~1の値で変わる.スコアを一般的な値で出力するならscale=0.負荷量のプロットをとスコアの相関係数の比率(値そのものは異なる)でプロットしたいならscale=1(これがデフォルト)

考え方

主成分分析のかみ砕いた解説(咀嚼説明)は,有馬・石村(1987)1にあります.

以下はそのフォローです.

成分ベクトルととの距離の最小化であり,

スコアの分散の最大化でもある.

情報損失量の最小化

と直線 との距離 を
で表す.

ただし, とする.

ヘッセの標準形から

ヘッセの標準形

から直線 に下した垂線の長さは

情報損失量 を以下で定義して,

を最小化するような を見つけたい.

言い換えれば・・・

s.t.

なので,ラグランジュ乗数法より

 (1)

1階の条件として

 (2)

なので,の平均を,それぞれ,とすると,

ということであり,とすると,

なので,

 (3)

 (4)

さらに,

 (5)

(3)式から

 (6)

(4)式から

 (7)

と書き換えると,

 (8)

ということであり,というのは,行2列の行列の分散・共分散行列の固有ベクトル

は,の分散・共分散行列の固有値

・・・ということ.

さらに,(6)式にを掛けて,(7)式にを掛けて足し合わせると

なので,

 (9)

であり,これは,

のことなので,情報損失量ということ.つまり,

 (10)

スコアの最大化

主成分得点 (principal component score)
あるいは単にスコア (score) と呼ぶ.

行2列の行列でスコアを最大化するノルム1のを求める.

であるから,ラグランジュ乗数法で


(10)

1階の条件は


(11)


(12)

したがって,


(13)

であり,は, の分散・共分散行列の固有ベクトル

は,その固有値.

なので

(14)

でもある.

Rで

Rのpcrcom関数で

データ

わかりやすいデータを作成.

の2変数9個のデータ.

library(tidyverse)
x01<-c(0,0,
  2+0.5,1-1,
  2-0.5,1+1,
  4+0.5,2-1,
  4-0.5,2+1,
  -2+0.5,-1-1,
  -2-0.5,-1+1,
  -4+0.5,-2-1,
  -4-0.5,-2+1
  )%>%
  matrix(ncol=2,byrow=T)

こんな感じ

plot(x01,asp=1)
grid()
text(1:9,x=x01[,1],y=x01[,2],pos = 1)
abline(h=0)
abline(v=0)
abline(a=0,b=0.5,col=4)
abline(a=0,b=-2,col=3)

横軸がx01[,1]),縦軸がx01[,2]

8つの点が,の方向の原点を通る直線対して垂直に,長さで対称にばらついている.

分析しなくても主成分は,の方向との方向であることはわかる.

平均は2変数とも0.

apply(x01,2,mean)%>%
  round(4)
## [1] 0 0

分散・共分散行列

var(x01)
##       [,1] [,2]
## [1,] 10.25  4.5
## [2,]  4.50  3.5

主成分分析を行う

prcomp()関数

pr01<-prcomp(x01)
pr01
## Standard deviations (1, .., p=2):
## [1] 3.535534 1.118034
## 
## Rotation (n x k) = (2 x 2):
##             PC1        PC2
## [1,] -0.8944272 -0.4472136
## [2,] -0.4472136  0.8944272

固有値

Standard deviationsというのは,スコアの標準偏差

pr01$sdev
## [1] 3.535534 1.118034

二乗するとスコアの分散であり,の分散・共分散行列の固有値

pr01$sdev^2
## [1] 12.50  1.25

この場合固有値は2個あるので,2つ出力.最初が第1主成分の,次のが第2主成分の固有値.第1主成分の方に分散が大きいことがわかる.

主成分のベクトル

Rotationは,主成分のベクトル

pr01$rotation
##             PC1        PC2
## [1,] -0.8944272 -0.4472136
## [2,] -0.4472136  0.8944272

これも2種類あって,1列目が第1主成分で,2列目が第2主成分で,直交している.

第1主成分がの方向で,第2主成分がの方向で,ノルムが1となっている.

pr01$rotation[,1]^2
## [1] 0.8 0.2

(横軸が,縦軸がの場合)左下方向と左上方向で,思っていたのと逆だが,どっち方向にスコアをとるかだけの問題で,本質的に違いはない.

各点のスコア

各点のスコアは,$xで取り出せる.

pr01$x
##             PC1       PC2
##  [1,]  0.000000  0.000000
##  [2,] -2.236068 -1.118034
##  [3,] -2.236068  1.118034
##  [4,] -4.472136 -1.118034
##  [5,] -4.472136  1.118034
##  [6,]  2.236068 -1.118034
##  [7,]  2.236068  1.118034
##  [8,]  4.472136 -1.118034
##  [9,]  4.472136  1.118034

スコアはのことだから,x01Rotationを掛け合わせても得られる.

x01%*%pr01$rotation
##             PC1       PC2
##  [1,]  0.000000  0.000000
##  [2,] -2.236068 -1.118034
##  [3,] -2.236068  1.118034
##  [4,] -4.472136 -1.118034
##  [5,] -4.472136  1.118034
##  [6,]  2.236068 -1.118034
##  [7,]  2.236068  1.118034
##  [8,]  4.472136 -1.118034
##  [9,]  4.472136  1.118034

各点のスコアのグラフ

第1主成分が横軸,第2主成分が縦軸

plot(pr01$x,asp=1)
grid()
abline(h=0,col=4)
abline(v=0,col=3)
text(1:9,x=pr01$x[,1],y=pr01$x[,2],pos=1)

スコアのプロットをマイナス45度回して,さらに+-をひっくり返した形.
(青い線がもとの真ん中を通る直線)

例えば,2番目の点は,傾き1/2の直線(第1主成分方向)に垂線を下した点がだから,第1主成分が

-sqrt(2^2+1^2)
## [1] -2.236068

第2主成分は,の点からだけずらした分だから,

-sqrt(0.5^2+(-1)^2)
## [1] -1.118034

スコアの分散は固有値

var(pr01$x)%>%
  round(4)
##      PC1  PC2
## PC1 12.5 0.00
## PC2  0.0 1.25

主成分負荷量

個々の点の特徴を,2つの主成分に分けて見れることはわかった.

それではという変数は,主成分から見てどんな変数か,はどうか,というのが主成分負荷量
(principal component loading).

因子負荷量 (factor
loading)と呼ぶ人もいるが,因子分析の用語で気持ち悪い.

単に負荷量 (loading) と呼ぶ場合も.

一般的な定義は

vinv01<-1/apply(x01,2,sd)
diag(vinv01)%*%pr01$rotation%*%diag(pr01$sdev)%>%
  round(4)
##         [,1]    [,2]
## [1,] -0.9877 -0.1562
## [2,] -0.8452  0.5345

これはとスコアとの相関係数.

3項目への変換は

からきている.

ふつうに相関係数を求めると.

cor(x01,pr01$x)%>%
  round(4)
##          PC1     PC2
## [1,] -0.9877 -0.1562
## [2,] -0.8452  0.5345

定義した負荷量と一致する.

なので,主成分負荷量はの値をとる.

変数は,第1主成分と-0.9877の相関があるが,第2主成分とは-0.1562しかない.

変数は,第1主成分と-0.8452の相関があり,第2主成分とも0.5345の相関がある.

biplot

主成分分析のスコアと主成分負荷量は,biplot()関数で出力できる.

biplot(pr01,asp=1,scale=1)%>%
grid()

スコアと負荷量は全く別の数値だが,無理やり1つのグラフに入れてある.
下と左の目盛がスコア用で,上と右が負荷量のための目盛.

1~9の数字が各点のスコアで,赤い矢印が負荷量・・・だと思ったが,何か数値が違う.

この図を理解するには,の特異値分解をすればよくわかる.

特異値分解

の特異値分解は,Rを使えばsvd()関数で簡単にできる.

svd01<-svd(x01)
svd01%>%
  lapply(.,function(x) round(x,4))
## $d
## [1] 10.0000  3.1623
## 
## $u
##          [,1]    [,2]
##  [1,]  0.0000  0.0000
##  [2,] -0.2236 -0.3536
##  [3,] -0.2236  0.3536
##  [4,] -0.4472 -0.3536
##  [5,] -0.4472  0.3536
##  [6,]  0.2236 -0.3536
##  [7,]  0.2236  0.3536
##  [8,]  0.4472 -0.3536
##  [9,]  0.4472  0.3536
## 
## $v
##         [,1]    [,2]
## [1,] -0.8944 -0.4472
## [2,] -0.4472  0.8944
u01<-svd01$u
d01<-svd01$d
v01<-svd01$v

もノルム1の直交するベクトル

norm(u01,"2")
## [1] 1
norm(v01,"2")
## [1] 1
var(u01)%>%
  round(4)
##       [,1]  [,2]
## [1,] 0.125 0.000
## [2,] 0.000 0.125

で各点のスコアとなる. ただし,

u01%*%diag(d01)%>%
  round(4)
##          [,1]   [,2]
##  [1,]  0.0000  0.000
##  [2,] -2.2361 -1.118
##  [3,] -2.2361  1.118
##  [4,] -4.4721 -1.118
##  [5,] -4.4721  1.118
##  [6,]  2.2361 -1.118
##  [7,]  2.2361  1.118
##  [8,]  4.4721 -1.118
##  [9,]  4.4721  1.118

・・・であるはずなのだが,biplot()関数では,スコアを以下で表し,

主成分負荷量を以下で表す.


ただし,は,の実数で,通常のスコアの定義は,の場合である.

一方,biplot()のスコアはで出しているようだ.

このの値は,biplot()関数の中で,引数scale=で与えることができる.

の場合

スコアは

u01%*%diag(d01^(1-0))%>%
  round(4)
##          [,1]   [,2]
##  [1,]  0.0000  0.000
##  [2,] -2.2361 -1.118
##  [3,] -2.2361  1.118
##  [4,] -4.4721 -1.118
##  [5,] -4.4721  1.118
##  [6,]  2.2361 -1.118
##  [7,]  2.2361  1.118
##  [8,]  4.4721 -1.118
##  [9,]  4.4721  1.118

・・・で,通常のスコアと一致する.

負荷量は

v01%*%diag(d01^0)%>%
  round(4)
##         [,1]    [,2]
## [1,] -0.8944 -0.4472
## [2,] -0.4472  0.8944

で,これは主成分ベクトルそのもの.

biplotで確認

biplot(pr01,asp=1,scale=0)%>%
grid()


スコアの値は正確.

負荷量は,定義した通りではあるのだが,90度で直交して,負荷量が相関係数を表していない.(だって,主成分ベクトルだから・・・)

は,第1主成分にくっついているので,と第1主成分は相関が高い!」・・・とか,負荷量を赤い矢印の傾きで判断したいのなら,相関係数の傾きと同じになってほしい.

相関関係の傾きは

c(-0.1562/-0.9877,0.5345/-0.8452)
## [1]  0.1581452 -0.6323947

の傾きは

c(-0.4472/-0.8944,0.8944/-0.4472)
## [1]  0.5 -2.0

・・・と,全然違う.(だって主成分ベクトルだから)

の場合

スコアは

u01%*%diag(d01^(1-1))%>%
  round(4)
##          [,1]    [,2]
##  [1,]  0.0000  0.0000
##  [2,] -0.2236 -0.3536
##  [3,] -0.2236  0.3536
##  [4,] -0.4472 -0.3536
##  [5,] -0.4472  0.3536
##  [6,]  0.2236 -0.3536
##  [7,]  0.2236  0.3536
##  [8,]  0.4472 -0.3536
##  [9,]  0.4472  0.3536

負荷量は,

v01%*%diag(d01^1)%>%
  round(4)
##         [,1]    [,2]
## [1,] -8.9443 -1.4142
## [2,] -4.4721  2.8284

biplotで確認すると

biplot(pr01,asp=1,scale=1)%>%
grid()

負荷量の傾きは・・・

c(-1.4142/-8.9443,2.8284/-4.4721)
## [1]  0.1581119 -0.6324546

相関係数の傾きとほぼ一緒.(値は違うけど)

の各要素の平方根をとったものをと表すと・・・

スコアは

u01%*%diag(d01^(1-0.5))%>%
  round(4)
##          [,1]    [,2]
##  [1,]  0.0000  0.0000
##  [2,] -0.7071 -0.6287
##  [3,] -0.7071  0.6287
##  [4,] -1.4142 -0.6287
##  [5,] -1.4142  0.6287
##  [6,]  0.7071 -0.6287
##  [7,]  0.7071  0.6287
##  [8,]  1.4142 -0.6287
##  [9,]  1.4142  0.6287

負荷量は

v01%*%diag(d01^0.5)%>%
  round(4)
##         [,1]    [,2]
## [1,] -2.8284 -0.7953
## [2,] -1.4142  1.5905

確認する

biplot(pr01,asp=1,scale=0.5)%>%
grid()

負荷量の傾きは

c(-0.7953/-2.8284,1.5905/-1.4142)
## [1]  0.2811837 -1.1246641

・・・微妙.

結局,負荷量を相関係数として理解するなら,でbiplotを描くのが最もそれに近い.ただし,目盛は違うので要注意・・・ということ.

寄与率

スコアの分散が固有値だが,分散の和における個々の分散(固有値)の割合は寄与率と呼ばれる.

固有値が第2主成分まで,2つある場合
第1主成分の寄与率は

2変数しかない場合,の分散()と一致するので,この寄与率は第1主成分がの分散を説明する割合ということになる.

スコアの分散は

var(pr01$x)%>%
  round(4)
##      PC1  PC2
## PC1 12.5 0.00
## PC2  0.0 1.25

固有値は

pr01$sdev^2
## [1] 12.50  1.25

寄与率は

pr01$sdev^2/sum(pr01$sdev^2)%>%
  round(4)
## [1] 0.90909091 0.09090909

第1主成分の寄与率が90.9%で,第2主成分の寄与率が9.09%ということ.

寄与率は,summary()関数で出力できる.

summary(pr01)
## Importance of components:
##                           PC1     PC2
## Standard deviation     3.5355 1.11803
## Proportion of Variance 0.9091 0.09091
## Cumulative Proportion  0.9091 1.00000

Proportion of Varianceが寄与率.それを大きい方から順に足していったのが累積寄与率.Cumulative Proportionがそれ.この場合2変数なので,主成分は高々2個.なので,第2主成分で累積寄与率は1になる.

標準化

主成分分析は,の各変数の単位をどうとるかで,出力される値が変わってしまうので,(よほど変動を変えたくない場合を除き)通常は全ての変数の分散を1に標準化してから主成分分析を行うことが一般的なようだ.

を標準化する.

scale()

x02<-scale(x01)
plot(x02,asp=1)
grid()
abline(h=0)
abline(v=0)
abline(a=0,b=1,col=4)
abline(a=0,b=-1,col=3)
text(1:9,x=x02[,1],y=x02[,2],pos=4)

平均

apply(x02,2,mean)
## [1] 0 0

分散

var(x02)
##           [,1]      [,2]
## [1,] 1.0000000 0.7513055
## [2,] 0.7513055 1.0000000

主成分分析

pr02<-prcomp(x02)
pr02
## Standard deviations (1, .., p=2):
## [1] 1.3233690 0.4986928
## 
## Rotation (n x k) = (2 x 2):
##             PC1        PC2
## [1,] -0.7071068 -0.7071068
## [2,] -0.7071068  0.7071068

2変数しかない場合,主成分のベクトルは,2本の45度線.

事前にを標準化しなくてもprcomp()関数に引数scale=Tを与えても同じ結果となる.

prcomp(x01,scale=T)
## Standard deviations (1, .., p=2):
## [1] 1.3233690 0.4986928
## 
## Rotation (n x k) = (2 x 2):
##             PC1        PC2
## [1,] -0.7071068 -0.7071068
## [2,] -0.7071068  0.7071068

スコアのプロット

plot(pr02$x,asp=1)
grid()
abline(h=0,col=4)
abline(v=0,col=3)
text(1:9,x=pr02$x[,1],y=pr02$x[,2],pos=1)

主成分負荷量

cor(x02,pr02$x)
##             PC1        PC2
## [1,] -0.9357632 -0.3526291
## [2,] -0.9357632  0.3526291

biplot

の場合

biplot(pr02,scale = 0)
grid()


負荷量(赤線)の値が,なんかちょっと違う気がするが・・・この辺りは謎.

の場合

biplot(pr02,scale = 1)
grid()

biplot(pr02,scale = 0.5)
grid()

寄与率

summary(pr02)
## Importance of components:
##                           PC1    PC2
## Standard deviation     1.3234 0.4987
## Proportion of Variance 0.8757 0.1244
## Cumulative Proportion  0.8757 1.0000

第1主成分の寄与率が,標準化していない場合よりも小さくなった.

主成分分析の結果は変数のスケールに応じて変わるということ.つまり同じデータでも,例えばキログラムで測るのか,それともトンで測るのかでも結果が変わるということ.

それはまずい,ということで,主成分分析では各変数の分散を1に揃えるのが一般的.

ただし,合計に意味があるような場合.例えば,入試の成績などでは,各科目の分散を1に揃えない方が,分析の意味づけがやりやすくなるはず,たぶん.


1 有馬哲・石村貞夫『多変量解析のはなし』東京図書,1987