目次

・多次元データを要約する意味

・「ばらつきの量」は「情報の量」

・多次元データの要約方法：主成分分析の原理

・回転後の座標：主成分得点

・元の軸と新しい軸の相関：主成分負荷量

・構成する分散の割合：寄与率

・Rで行う主成分分析

・実データの主成分分析

・追記1 : ラグランジュの未定乗数法を用いた解説

・追記2 : 主成分得点から元のデータを復元する

多次元データを要約する意味

　データ解析では、極めて多数の変数を扱うこともあるだろう。特に、野外調査なんかでは環境データがその典型例で、ある地点の変数には、気温、湿度、酸素量、二酸化炭素量、日射量、土壌水分含量、窒素含量、地温……などと、挙げていけばきりがない。様々な地点の植物の成長量に影響量を与える変数を探索するために、これらの膨大な量の説明変数を盛り込むことは現実的ではない。統計モデリングを実行するためには、推定したいパラメータ数よりも被説明変数のサンプルサイズが大きい必要がある。けれども、こういう多数の変数を含む環境データをすべて説明変数に組み込むと、植物の成長量のデータに対して説明変数の数が膨大すぎて、解析は失敗に終わる。そもそも、これらの説明変数はすべて入力する必要はあるのだろうか。例えば、気温と地温は極めて高い相関を示すだろう。日射量と気温もそうだし、湿度と気温だって相関するだろう。湿度と土壌水分含量も関係がありそうだ。酸素量や二酸化炭素量なんて各環境で大きな違いはないだろう……となったとき、強い相関を示す変数はどれかを代表して除去するべきだし、ばらつきがあまりに小さい変数を解析に組み込む意味はない。あるいは、そもそもの目的として植物の成長量と関係のある変数を選びたいのだから、植物の成長量と相関の小さい変数を解析に組み込む意味もないだろう。こういう風に、多次元データをそのまま使うのではなく、数を絞る＝要約することで、データ解析もしやすくなるし、解釈もしやすくなる。

「ばらつきの量」は「情報の量」

　本稿では多次元データの要約の方法を紹介するが、その前に基本的なコンセプトを身に着けてもらいたい。例えば、下記のような2次元データがあったとする。この時、X1軸かX2軸のどちらかを切り捨ててデータの情報を抽出せよ＝要約せよ、といわれたらどうするだろう。データの情報っていうのは結構、抽象的な概念だが、イメージをつかむために考えてみよう。

ここで、X1軸を切り捨ててX2軸を選択することを考える。すると全くばらつきのないデータが得られるだろう。ばらつきがない、というのも情報だけど、もはやその情報しか残っていない。本当は切り捨てたX1軸にはばらついていたのに、X2軸でデータを抽出したためにばらつきの情報がすべて消えて、あたかも一定の値になってしまっている。逆にX2を切り捨ててもデータには何の影響なく、X1軸に沿った情報は元のデータと完全に一致している。つまり、X1軸の情報で元のデータを完全に復元できる。このように、ばらつきが大きい軸というのは多くの情報を含んでいる、と解釈が可能である。

多次元データの要約方法：主成分分析の原理

　データを要約するうえでは、元の情報をなるべく維持しつつも、重要ではない変数をそぎ落とす必要がある。そこで、多次元データをプロットした時に、上記のコンセプトにならって、最もばらついている軸を探索することを考える。これは元の軸からばらつきの大きい軸を選ぶということではなく、新たに最もばらついている方向に軸を設定するという意味である。この方法を主成分分析Principal component analysis（PCA）と呼び、非常によく使われるデータの要約方法である。近年の機械学習の文脈では、教師なし学習の1種として紹介されることもある（一方で、線形回帰などは教師あり学習の1種。機械学習の分野で一見、小難しいことを言っているようで、ただの主成分分析でした……みたいなことが頻発しているように思う）。簡単のため2次元データに基づいた解釈を紹介していくが、多次元になっても同様の考え方でよい。

もっともデータがばらつく軸を探索するということなので、多次元上でのある軸に沿ったデータのばらつきを定義する必要がある。まず、任意の軸に対するデータのばらつきとは下図の青矢印の二乗和を指す。

なので、青矢印＝あるデータ点を任意の軸に射影した長さを算出する。下図のように計算ができる。

ゆえに、あるデータ点を任意の軸に射影した長さの二乗和は下図のようになり、それをサンプルサイズで割ったものを任意の軸方向の分散とする。その後、以下のように式変形できる。

以上の式変形から、任意の軸方向の分散は軸方向の単位ベクトルuと分散共分散行列Σを使って表すことができた。このような(ベクトル)^T(行列)(ベクトル)というような表示形式を二次形式quadratic formと呼ぶことがある。

　では、この分散が最大になるようなuを求める。実はこの問題は分散共分散行列Σの固有ベクトルを求める問題に帰着ができる。本当は、ある定理を認めれば比較的早く証明できるが（追記にて紹介）、ここでは純朴に式変形のみで示す。激しい式変形が続くが、やっていることは関数の最大値を求める問題にすぎない。

uとなる具体的なφやsin(φ/2)、cos(φ/2)は求められないが、軸の傾きを求めるだけならsinφ、cosφから計算ができる。分散を最大にする軸の傾きは(-cosφ+1)/sinφ、最小にする傾きは(-cosφ-1)/sinφとなる。また、あわせてそのときの分散の最大値・最小値=f(θ)も覚えておいてほしい。

　一方で、分散共分散行列の固有ベクトルeigenvectorは以下のようにして求めることができる。ここでは固有ベクトルの求め方を解説しないので、必要に応じて復習してほしい。計算は重たいが、最終的にはシンプルな形に落ち着く。

このように2つの異なる方法で求めたu（軸の傾き）が一致する。なお、固有ベクトルとして求めると分散を最大にする軸と最小にする軸が同時に求まるが、最大にする軸は絶対値が大きい方の固有値eigenvalue(λ)に対応している。ここで固有値λを見ると、これは軸方向の分散f(θ)と一致している。以上の考察から、データの分散が最大となるような軸とその時の分散は、分散共分散行列の固有値・固有ベクトルを求めればよいことがわかる。また、求めた2つの軸は直交する。これは(最大にする傾き)×(最小にする傾き)=-1や(軸の方向のベクトルの内積)=0となるので、すぐに確認できる。次元が増えても「最大の固有値と対応する固有ベクトルが分散を最大にする軸」「固有値は対応する固有ベクトル方向の分散」「固有ベクトル（求まった軸）同士は直交する」という関係は変わらない。PCAによって重要な＝分散の大きい方向の軸を選抜する場合、固有値が大きいものから選択すればよい。固有ベクトルが直交する性質から、後述する主成分得点を説明変数とすることで、多重共線性の問題をクリアすることもできる。

　実際のデータ解析では、分散共分散行列をそのまま使うことは少ない。データの分散は当然、データの値のオーダーに影響を受けるためである。そこで、データの値を平均0、分散1にスケーリングしてから分散共分散行列を作り、それの固有値・固有ベクトルを求める。このようなスケーリング後の分散共分散行列は、x1とx2の相関係数r12を用いて以下のようにあらわせ、これを相関行列correlation matrixと呼ぶ。

つまり、主成分分析は多くの場合、相関行列の固有値・固有ベクトルを求めていることになる。上記で求まった条件にVar(x1) = 1、Var(x2) = 1、Cov(x1, x2) = r12として計算していけば以下のように求まる。

以上のように非常にシンプルな形で表せ、分散=固有値λ = 1 + r12 or 1 - r12となった。

回転後の座標：主成分得点

　ここまでで分散が最も大きい軸と小さい軸を固有ベクトルで表せることを示してきた。この固有ベクトルを新たな軸として図を描きなおすことを考える。下図のように固有ベクトルの傾き分だけ逆回転すれば、固有ベクトルを軸とした図に描き換えることができる。

このような変換後の座標を主成分得点principal component scoreと呼ぶ。また、この時新たな軸となる固有ベクトルの方向を、固有値の絶対値の大きい方から、第1主成分（principal component 1、PC1）、第2主成分（PC2）……と名前が付けられる。主成分得点は図の様子から、平均を引いた元のデータと固有単位ベクトルの内積をとれば計算が可能である。主成分分析時には、この主成分得点を描画する。逆に主成分得点と固有ベクトルの情報があれば、元のデータを復元することもできる（追記）。

元の軸と新しい軸の相関：主成分負荷量

　主成分分析のもともとの目的は変数の選択であった。固有ベクトルをもとに貼られた空間に対して、元の軸がどれくらい相関しているかを明らかにできれば、それをもとに変数の選択ができそうだ。このような元の軸方向のデータと新しい軸方向のデータの相関の強さを主成分負荷量principal component loading（あるいは因子負荷量とも）と呼ぶ。主成分負荷量は、下図の回転後の空間におけるベクトルで表されることが多い。例えば、X1の主成分負荷量はX1の値とPC1との相関、PC2との相関……をベクトルとして示したものである。下図のように、相関を直接計算しなくても固有値の平方根と固有ベクトルの要素の積で計算ができることが知られている。

下図のように、主成分負荷量をもとの変数の標準偏差で除した値が、相関係数となる。この相関係数の形になっているものを主成分負荷量とする場合もある。また、解析した行列が相関行列の場合、変数の標準偏差を1に統一しているので、主成分負荷量の値は相関係数そのものになる。

主成分負荷量を示すベクトルが相対的に長いほど、その変数はその方向の主成分との相関が高いことを示す。つまりは、主成分負荷量の長さが長い変数を選択すれば、より効率よくデータを要約できる、ということである。また、各変数の主成分負荷量の示すベクトルの方向が大体同じであれば、その変数間では正の相関があり、逆に主成分負荷量の方向が逆であれば、その変数間では負の相関があることがわかる。もし、主成分負荷量のベクトルが一直線に並ぶ場合は、その変数間には極めて強い相関があることが疑われる。

構成する分散の割合：寄与率

　一方で、主成分の分散はデータ全体の分散のどれだけの割合を占めるのか。この割合のことを寄与率contribution ratioと呼ぶ。寄与率が大きければ大きいほど、その主成分はデータのばらつきの多くの原因となっていることがわかる。この計算はシンプルで、(ある主成分の固有値)/(固有値＝分散の総和)で表すことができる。通常は、下図のように主成分とともに寄与率を記載する。下図では、第1主成分が89.8%の分散を占め、第2主成分が10.2%を占める。

Rで行う主成分分析

　長々と解説してきたが、Rでの実際の解析はそれほど難しいものではない。ついでに、いくつか手計算して、値が一致することを確認しよう。以下のようにデータを生成し、その様子を作図してみる。

------------------------------------------------------

library(plotn)

library(MASS)

s1 <- 20

s2 <- 5

s12 <- 5

n <- 10000

mu <- c(0,0)

v <- matrix(c(s1, s12, s12, s2), ncol = 2)

v

##      [,1] [,2]

## [1,]   20    5

## [2,]    5    5

d <- as.data.frame(mvrnorm(n, mu = mu, Sigma = v)) #相関するデータの生成

colnames(d) <- c("X1","X2")

head(d)

##           X1         X2

## 1  5.9839318 -0.2709100

## 2  2.6880302  1.0835966

## 3 -4.8870678 -1.4494175

## 4 -5.3332774 -0.1232362

## 5 -0.7712681  3.3682141

## 6 -3.3904722 -0.6370487

plotn(X2 ~ X1, data = d, col.dot = "#00000020", asp = 1) #図2の作図

------------------------------------------------------

では以下に、主成分分析を実行してみる。prcomp関数に多変量データをインプットすれば、解析結果を出してくれる。今回は、スケーリングをしないで解析してみよう。以下に一通りの解析結果の表示方法を載せる。寄与率に関しては、Proportion of Varianceの値が寄与率である。

------------------------------------------------------

#以下、主成分分析のコマンド

res <- prcomp(x = d, scale = F) #主成分分析、スケーリングをしない

res$sdev #固有値の平方根、標準偏差

## [1] 4.661726 1.897214

(res$sdev)^2 #固有値、分散

## [1] 21.731693  3.599422

res$rotation #固有ベクトル

##           PC1        PC2

## X1 -0.9556585  0.2944773

## X2 -0.2944773 -0.9556585

head(res$x) #主成分得点

##             PC1          PC2

## [1,] -5.6981069  2.065900570

## [2,] -2.9472221 -0.199113420

## [3,]  5.0378997 -0.009111671

## [4,]  5.0737933 -1.407886587

## [5,] -0.3140824 -3.401112426

## [6,]  3.3684412 -0.344745290

b1 <- res$rotation[2,1]/res$rotation[1,1] #第1主成分の傾き

b2 <- res$rotation[2,2]/res$rotation[1,2] #第2主成分の傾き

b1

## [1] 0.3081408

b2

## [1] -3.24527

l <- t(t(res$rotation)*res$sdev) #主成分負荷量

l

##          PC1        PC2

## X1 -4.455018  0.5586866

## X2 -1.372773 -1.8130890

l2 <- l*c(1/sd(d$X1), 1/sd(d$X2)) #主成分負荷量（相関係数）

l2

##           PC1        PC2

## X1 -0.9922282  0.1244315

## X2 -0.6036399 -0.7972571

lo <- summary(res) #寄与率

lo

## Importance of components:

##                           PC1    PC2

## Standard deviation     4.6617 1.8972

## Proportion of Variance 0.8579 0.1421

## Cumulative Proportion  0.8579 1.0000

------------------------------------------------------

では、これらをすべて手計算してみる。以下のように計算できる。eigen関数は行列の固有値・固有ベクトルを解く関数である。値が一通り、一致していることを確認しよう。ただし、たまに固有ベクトルが逆方向にとられることがあり、主成分得点が異符号になることがある。しかし、実質的な違いはない。

------------------------------------------------------

#以下、手計算での主成分分析

cm <- cov(d) #分散共分散行列

cm

##           X1       X2

## X1 20.159318 5.102779

## X2  5.102779 5.171797

e <- eigen(cm) #固有値・固有ベクトルを解く

e$values #固有値、分散

## [1] 21.731693  3.599422

e$vectors #固有ベクトル

##            [,1]       [,2]

## [1,] -0.9556585  0.2944773

## [2,] -0.2944773 -0.9556585

dt <- d

dt$X1 <- d$X1 - mean(d$X1) #平均値だけずらす

dt$X2 <- d$X2 - mean(d$X2) #平均値だけずらす

x <- as.matrix(dt)%*%e$vectors #主成分得点

colnames(x) <- c("PC1","PC2")

head(x)

##             PC1          PC2

## [1,] -5.6981069  2.065900570

## [2,] -2.9472221 -0.199113420

## [3,]  5.0378997 -0.009111671

## [4,]  5.0737933 -1.407886587

## [5,] -0.3140824 -3.401112426

## [6,]  3.3684412 -0.344745290

s1e <- var(d$X1)

s2e <- var(d$X2)

s12e <- cov(d$X1,d$X2)

sin_e <- (2*s12e)/sqrt((s1e-s2e)^2+4*s12e^2) #sin

cos_e <- (s1e-s2e)/sqrt((s1e-s2e)^2+4*s12e^2) #cos

slope <- (-cos_e + c(1,-1))/sin_e #各主成分の傾き

slope

## [1]  0.3081408 -3.2452702

cor(d, x) #主成分負荷量（相関係数）

##           PC1        PC2

## X1 -0.9922282  0.1244315

## X2 -0.6036399 -0.7972571

e$values/sum(e$values) #寄与率

## [1] 0.8579051 0.1420949

------------------------------------------------------

では、上記の解析結果を図示しよう。まず、元の図に固有ベクトルを描きこんでみよう。

------------------------------------------------------

plotn(X2 ~ X1, data = d, col.dot = "#00000020", asp = 1, xlab = "X1", ylab = "X2") #図3の描画

overdraw(arrows(mean(d$X1), mean(d$X2), 

                res$sdev[1]*res$rotation[1,1]+mean(d$X1),

                res$sdev[1]*res$rotation[2,1]+mean(d$X2), 

                col = "blue", length = 0.1),

         arrows(mean(d$X1), mean(d$X2),

                res$sdev[2]*res$rotation[1,2]+mean(d$X1), 

                res$sdev[2]*res$rotation[2,2]+mean(d$X2), 

                col = "blue", length = 0.1))

------------------------------------------------------

青矢印で固有ベクトルを示し、長さは標準偏差に対応している。より長い方の固有ベクトルが分散が大きい方向に向いていることを確認できる。また、固有ベクトルは互いに直交している。続いて、主成分得点と主成分負荷量を図示しよう。この図示が、主成分分析の結果としてよく示される。なお、以下のように複雑なコマンドを使ったが、簡単に確認するだけなら後述するbiplot関数もある（きれいな図ではないので、いずれにしても自分なりに作図が必要になる）。

------------------------------------------------------

new_d <- as.data.frame(res$x)

plotn(PC2 ~ PC1, data = new_d, col.dot = "#00000020", asp = 1, 

      xlab = paste0("PC1 (", round(lo$importance[2,1]*100,digits = 1),"%)"), 

      ylab = paste0("PC2 (", round(lo$importance[2,2]*100,digits = 1),"%)")) #図4の描画

overdraw(arrows(0, 0, l[1,1], l[1,2], length=0.1, col="red"),

         arrows(0, 0, l[2,1], l[2,2], length=0.1, col="red"))

------------------------------------------------------

図3と図4を比べると、図3の長い固有ベクトルを右方向に向くまで回転した図が図4になっていることがわかる。寄与率を確認すると、第1主成分が分散の85.8%を占め、第2主成分が14.2%を占めている。ずまり、ばらつきのほとんどは第1主成分を見るだけで理解できる、ということである。赤矢印で主成分負荷量を示している。第1主成分とX1の主成分負荷量がほぼ反対の方向に向き、長いので第1主成分とX1は高い逆相関の関係になっていることが読み取れる。事実、図3を確認すると、長い方の固有ベクトルとX1が逆向きであり、固有ベクトルもほぼ水平なので高い相関を示すことがわかる。一方、X2は主成分負荷量が短く、第1主成分とほぼ垂直方向にあり、相関が高くないことが読み取れる。これも、図3から読み取れることと同様である。また、X1とX2の主成分負荷量はほぼ同じ方向を向いている（独立に近いが）ので正の相関があることがわかるが、同様にこれも図3から読み取れる。

　今度はスケーリングした場合の解析例を示そう。以下のようにデータを生成する。今度は負の相関があることに注目しよう。

------------------------------------------------------

s1 <- 20

s2 <- 5

s12 <- -5

n <- 10000

mu <- c(0,0)

v <- matrix(c(s1, s12, s12, s2), ncol = 2)

v

##      [,1] [,2]

## [1,]   20   -5

## [2,]   -5    5

d <- as.data.frame(mvrnorm(n, mu = mu, Sigma = v)) #相関するデータの生成

colnames(d) <- c("X1","X2")

ds <- data.frame(Scale_X1 = scale(d$X1), Scale_X2 = scale(d$X2)) #データのスケーリング

head(ds)

##     Scale_X1     Scale_X2

## 1 -2.1036738 1.4363113360

## 2  0.9309951 0.3844354528

## 3 -1.1544973 0.0647787799

## 4  0.7648902 0.1182638320

## 5 -0.3772391 0.3969752629

## 6  0.6805838 0.0009677197

plotn(Scale_X2 ~ Scale_X1, data = ds, col.dot = "#00000020", asp = 1, 

      xlab = "Scaled X1", ylab = "Scaled X2") #図5の作図

------------------------------------------------------

では、主成分分析を行おう。prcomp関数はscale = Tとすれば、自動でスケーリングしてくれるので、スケーリング前のデータセットを入力してよい。

------------------------------------------------------

#以下、主成分分析のコマンド

res <- prcomp(x = d, scale = T) #主成分分析、スケーリングする

res$sdev #固有値の平方根、標準偏差

## [1] 1.2230688 0.7100019

(res$sdev)^2 #固有値、分散

## [1] 1.4958973 0.5041027

res$rotation #固有ベクトル

##           PC1        PC2

## X1  0.7071068 -0.7071068

## X2 -0.7071068 -0.7071068

head(res$x) #主成分得点

##             PC1         PC2

## [1,] -2.5031475  0.47189652

## [2,]  0.3864760 -0.93014987

## [3,] -0.8621584  0.77054738

## [4,]  0.4572339 -0.62448421

## [5,] -0.5474523 -0.01395555

## [6,]  0.4805611 -0.48192969

b1 <- res$rotation[2,1]/res$rotation[1,1] #第1主成分の傾き

b2 <- res$rotation[2,2]/res$rotation[1,2] #第2主成分の傾き

b1

## [1] -1

b2

## [1] 1

l <- t(t(res$rotation)*res$sdev) #主成分負荷量（相関係数）

l

##           PC1        PC2

## X1  0.8648402 -0.5020472

## X2 -0.8648402 -0.5020472

lo <- summary(res) #寄与率

lo

## Importance of components:

##                          PC1   PC2

## Standard deviation     1.223 0.710

## Proportion of Variance 0.748 0.252

## Cumulative Proportion  0.748 1.000

------------------------------------------------------

つづいて手計算で主成分分析してみる。

------------------------------------------------------

#以下、手計算での主成分分析

cm <- cor(d) #相関行列

cm

##            X1         X2

## X1  1.0000000 -0.4958973

## X2 -0.4958973  1.0000000

e <- eigen(cm) #固有値・固有ベクトルを解く

e$values #固有値、分散

## [1] 1.4958973 0.5041027

1+c(1,-1)*cor(d$X1,d$X2) #固有値、分散、相関係数を使って計算もできる

## [1] 0.5041027 1.4958973

e$vectors #固有ベクトル

##            [,1]       [,2]

## [1,] -0.7071068 -0.7071068

## [2,]  0.7071068 -0.7071068

x <- as.matrix(ds)%*%e$vectors #主成分得点

colnames(x) <- c("PC1","PC2")

head(x)

##             PC1         PC2

## [1,]  2.5031475  0.47189652

## [2,] -0.3864760 -0.93014987

## [3,]  0.8621584  0.77054738

## [4,] -0.4572339 -0.62448421

## [5,]  0.5474523 -0.01395555

## [6,] -0.4805611 -0.48192969

cor(ds, x) #主成分負荷量（相関係数）

##                 PC1        PC2

## Scale_X1 -0.8648402 -0.5020472

## Scale_X2  0.8648402 -0.5020472

e$values/sum(e$values) #寄与率

## [1] 0.7479486 0.2520514

------------------------------------------------------

ではこれを図示しよう。まずは、元データに固有ベクトルを描きこむ。

------------------------------------------------------

plotn(Scale_X2 ~ Scale_X1, data = ds, col.dot = "#00000020", asp = 1, 

      xlab = "Scaled X1", ylab = "Scaled X2") #図6の作図

overdraw(arrows(mean(d$X1), mean(d$X2), 

                res$sdev[1]*res$rotation[1,1]+mean(d$X1),

                res$sdev[1]*res$rotation[2,1]+mean(d$X2), 

                col = "blue", length = 0.1),

         arrows(mean(d$X1), mean(d$X2),

                res$sdev[2]*res$rotation[1,2]+mean(d$X1), 

                res$sdev[2]*res$rotation[2,2]+mean(d$X2), 

                col = "blue", length = 0.1))

------------------------------------------------------

では、主成分得点と主成分負荷量を図示しよう。

------------------------------------------------------

new_d <- as.data.frame(res$x)

plotn(PC2 ~ PC1, data = new_d, col.dot = "#00000020", asp = 1, 

      xlab = paste0("PC1 (", round(lo$importance[2,1]*100,digits = 1),"%)"), 

      ylab = paste0("PC2 (", round(lo$importance[2,2]*100,digits = 1),"%)")) #図7の作図

for(i in 1:nrow(l)){

  overdraw(arrows(0, 0, l[i,1], l[i,2], length=0.1, col="red"))

  if (l[i,1] >= 0) {

    angle <- atan2(l[i,2], l[i,1]) / pi * 180

    text(l[i,1], l[i,2], labels = row.names(l)[i], 

         adj = -0.1, col= "red", srt = angle)

  } else {

    angle <- atan2(l[i,2], l[i,1]) / pi * 180 + 180

    text(l[i,1], l[i,2], labels = row.names(l)[i], 

         adj = +1.1, col = "red", srt = angle)

  }

}

------------------------------------------------------

今度はX1とX2の主成分負荷量が反対方向を向いているので、X1とX2が負の相関を示すことがわかる。

実データの主成分分析

　とはいえ、実データがないと主成分分析のありがたみはわかりにくい。そこで、Rにプリセットされているデータを使って解析をしてみよう。data("iris")と入力するとirisというオブジェクトが作られる。この中にはアヤメ属の3種の萼片の長さと幅、花弁の長さと幅が入っている。このデータセットは極めて有名なもので、統計学で功績を残したRonald Aylmer Fisher (1936) が論文で使用したものとされる。なお、実際にデータ収集したのは、Edgar Anderson (1935)である。

　まず、このデータセットの萼片の長さと幅を使って作図してみよう。

------------------------------------------------------

data("iris")

?iris #irisデータセットのヘルプ

d <- iris

head(d)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species

## 1          5.1         3.5          1.4         0.2  setosa

## 2          4.9         3.0          1.4         0.2  setosa

## 3          4.7         3.2          1.3         0.2  setosa

## 4          4.6         3.1          1.5         0.2  setosa

## 5          5.0         3.6          1.4         0.2  setosa

## 6          5.4         3.9          1.7         0.4  setosa

plotn(Sepal.Width ~ Sepal.Length, data = d, mode = "m", group = "Species", 

      legend = T, leg.sp = 5) #図8の描画

------------------------------------------------------

見ると、setosaはほかの2種と異なる傾向があることがわかる。実際には4つの測定項目があり、描画できる組み合わせは6通り考えられる。興味があればほかのセットでも確認してみよう。これを一度に見て、データの性質を見極めるのは困難だろう。こんな時に、データの要約として主成分分析が役立つ。

------------------------------------------------------

res <- prcomp(x = d[,-5], scale = T)

res

## Standard deviations (1, .., p=4):

## [1] 1.7083611 0.9560494 0.3830886 0.1439265

## 

## Rotation (n x k) = (4 x 4):

##                     PC1         PC2        PC3        PC4

## Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863

## Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096

## Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492

## Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

res$sdev

## [1] 1.7083611 0.9560494 0.3830886 0.1439265

res$rotation

##                     PC1         PC2        PC3        PC4

## Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863

## Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096

## Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492

## Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

l <- t(t(res$rotation)*res$sdev) # 主成分負荷量

l

##                     PC1         PC2         PC3         PC4

## Sepal.Length  0.8901688 -0.36082989  0.27565767  0.03760602

## Sepal.Width  -0.4601427 -0.88271627 -0.09361987 -0.01777631

## Petal.Length  0.9915552 -0.02341519 -0.05444699 -0.11534978

## Petal.Width   0.9649790 -0.06399985 -0.24298265  0.07535950

lo <- summary(res) #寄与率

lo

## Importance of components:

##                           PC1    PC2     PC3     PC4

## Standard deviation     1.7084 0.9560 0.38309 0.14393

## Proportion of Variance 0.7296 0.2285 0.03669 0.00518

## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

------------------------------------------------------

これを作図してみよう。

------------------------------------------------------

new_d <- data.frame(as.data.frame(res$x), Species = d[,5])

plotn(PC2 ~ PC1, data = new_d, group = "Species", mode = "m", asp = 1, 

      xlab = paste0("PC1 (", round(lo$importance[2,1]*100,digits = 1),"%)"), 

      ylab = paste0("PC2 (", round(lo$importance[2,2]*100,digits = 1),"%)"), 

      legend = T, leg.sp = 5) #図9の描画

for(i in 1:nrow(l)){

  overdraw(arrows(0, 0, l[i,1], l[i,2], length=0.1, col="red"))

  if (l[i,1] >= 0) {

    angle <- atan2(l[i,2], l[i,1]) / pi * 180

    text(l[i,1], l[i,2], labels = row.names(l)[i], 

         adj = -0.1, col= "black", srt = angle)

  } else {

    angle <- atan2(l[i,2], l[i,1]) / pi * 180 + 180

    text(l[i,1], l[i,2], labels = row.names(l)[i], 

         adj = +1.1, col = "black", srt = angle)

  }

}

------------------------------------------------------

まず、第1主成分が73%を占めているので、データセットのばらつきはこの軸方向で73%占めていることがわかる。また、各種の分布をみると、第1主成分に沿って分かれている。なので、この軸に沿った主成分負荷量が特に重要な意味を持っていると判断できる。第1主成分とPetal.Length（花弁の長さ）・Petal.Width（花弁の幅）が並行で、これらの変数がばらつきと強く関係している。また、花弁の長さと花弁の幅は極めて強い正の相関がある。もし、これらの変数を説明変数にするなら、どちらかを選択すればよいだろう。また、Sepal.Length（萼片の長さ）も第1主成分と強く関係していることがわかる。これら3つの変数があれば、3種はおおむね区別可能になる、というのが主成分分析の結果となる。一方で、Sepal.Width（萼片の幅）は種の識別にあまり役立たないことがわかる。実際、図8に立ち返れば、それが理解できるだろう。なお、上のスクリプトはきれいな図を作るために煩雑なものになっているが、簡単に確認するだけなら、下記のようなbiplot関数でもよい。ただし、本当に確認用である。

------------------------------------------------------

biplot(res)

------------------------------------------------------

追記1 : ラグランジュの未定乗数法を用いた解説

　ある関数f(x1, x2……xm)について最大・最小値を求めることを考えるとき、この関数を偏微分して、偏微分係数∂f/∂x1 = 0、∂f/∂x2 = 0……∂f/∂xm = 0を満たすx1、x2……xmの組を求めることで極値を調べるのが一般的である（偏微分に関しては線形モデルの項を参照）。しかし、g(x1, x2……xm) = 0となるようにx1……xmがとる値に制約が課されている場合、関数fの偏微分係数を求めるだけでは最大・最小値を探すことはできない。そこで、等式制約条件gがある中で、ある関数fの最大・最小値を求めるときによく用いられる方法がラグランジュの未定乗数法method of Lagrange multiplierと呼ばれるものである。

　ラグランジュ乗数λという未知数を導入し、以下の関数について各変数x1……xmおよびλの偏微分を考える。

このとき、以下の等式を満たすようなx1、x2……xm、λの組が存在し、この組が等式制約条件gがある中での関数fの最大・最小値を与える候補点である、というものである。

この性質を利用し、ベクトルuに関して下記の制約条件g(u1,u2……um) = 0を課した時の、関数f(u1,u2……um)の最大値を与える候補点を考える。

ここで一度、ベクトルと行列の積で表される関数f(u1,u2……um)について、u1に注目しながら展開する。CおよびC1～Cm-1はu1を含まない定数としてある。

この関数fをu1で偏微分すると以下のようになる。

ところで関数fはu1～umに関して対称な式だから、同様に考えれば関数fのu1～umに関する偏微分は以下のようになる。

一方、制約条件の関数gに関するu1～umの偏微分は以下のようになる。

以上から、関数Fのu1～umの偏微分は以下のようになる。式変形を続けると、uが満たすべき条件にとして分散共分散行列の固有値・固有ベクトルを求める式に帰着する。

また、上記の式を変形することで固有値が固有ベクトル方向の分散をであることを示せる。

なぜ、ラグランジュの未定乗数法が妥当な方法なのかというと、g(x1, x2, ……)の制約条件の下、f(x1, x2……)が最大値をとるときは、その最大値の与える候補点（x1', x2'……）においてfとgの接線は共有されているはずである。以下の図のように、2変数で考えれば、接線が共有されているということは、接線に垂直なベクトルも共有されている。

接線に垂直なベクトルは簡単に求められて、元の関数を各変数で偏微分したものを要素として持っている。これを勾配ベクトルと呼び、ある点での関数の最も大きな変化を示す方向を指す。ここでの主張は2つの関数fとgが勾配ベクトルを共有するということだが、方向は等しくても2関数のベクトルの長さが同じとは限らない。そこで、λ倍だけ調整することを考えて、以下の等式を満たせばよく、式変形すればラグランジュの未定乗数法に帰着される。

追記2 : 主成分得点から元のデータを復元する

　元データxiから固有ベクトルeと主成分得点xi'が計算されるということは、逆に固有ベクトルと主成分得点があれば元のデータを復元できる、ということである。以下のように計算すればよいことがわかる（平均値だけは主成分分析時に消失しているが）。

コマンドでは以下のようにできる。

------------------------------------------------------

s1 <- 20

s2 <- 5

s12 <- 5

n <- 10000

mu <- c(0,0)

v <- matrix(c(s1, s12, s12, s2), ncol = 2)

d <- as.data.frame(mvrnorm(n, mu = mu, Sigma = v)) #相関するデータの生成

colnames(d) <- c("X1","X2")

res <- prcomp(x = d, scale = F) #主成分分析、スケーリングをしない

res$rotation #固有ベクトル

##           PC1        PC2

## X1 -0.9556585  0.2944773

## X2 -0.2944773 -0.9556585

head(res$x) #主成分得点

##             PC1          PC2

## [1,] -5.6981069  2.065900570

## [2,] -2.9472221 -0.199113420

## [3,]  5.0378997 -0.009111671

## [4,]  5.0737933 -1.407886587

## [5,] -0.3140824 -3.401112426

## [6,]  3.3684412 -0.344745290

#逆算で元のデータを復元

head(res$x %*% solve(res$rotation) + rep(colMeans(d), each = nrow(d)))

##           X1         X2

## 1  5.9839318 -0.2709100

## 2  2.6880302  1.0835966

## 3 -4.8870678 -1.4494175

## 4 -5.3332774 -0.1232362

## 5 -0.7712681  3.3682141

## 6 -3.3904722 -0.6370487

head(d)

##           X1         X2

## 1  5.9839318 -0.2709100

## 2  2.6880302  1.0835966

## 3 -4.8870678 -1.4494175

## 4 -5.3332774 -0.1232362

## 5 -0.7712681  3.3682141

## 6 -3.3904722 -0.6370487

------------------------------------------------------

挑戦的なことを言えば、固有値の絶対値が小さい固有ベクトルを除去することで、近似的なデータを逆生成することができる。