Rによる主成分分析
主成分分析 (principal component analysis) とは多次元のデータを低次元データに縮約する方法のことである．PCA とも呼ばれる．高次元データを2次元か3次元に落とすことで人間が理解しやすい形式に変換するために行う．R で主成分分析を行う関数には princomp() と prcomp() の2種類が存在するが，princomp() にはサンプルサイズと変数の数に制限があるので，prcomp() を使えば，その点を意識せずに済む．以下に示されるような R に付属の11次元からなるデータ mtcars に対して主成分分析を行い次元の縮約を行う．mtcars の次元は，燃費 mpg，シリンダーの数 cyl，排気量 displacement，馬力 hp，重量 wt，ドラッグレースのタイム qsec，V型エンジンか直列型エンジンかの分類 vs，トランスミッション am，ギア数 gear，キャブレーター数 carb の11個で，それに対して，26個の車のデータが格納されている．
mpg cyl disp hp drat wt qsec vs am gear carb Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4 Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4 Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1 Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1 Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2 Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1 Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4 Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2 Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2 Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4 Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4 Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3 Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3 Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3 Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4 Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4 Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4 Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1 Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2 Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1 Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1 Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2 AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2 Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4 Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2 Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
このデータは以下のコマンドで確認できる．
$ mtcars
主成分分析の実行
これに対して PCA を行うには以下のように打つ．これによりオブジェクト rpca に主成分分析の結果が格納される．ここで，scale=T は数値のスケールが合っていないときに用いるオプションであり，スケールが統一されている場合は入力する必要はない．スケールが合っていないとは，次元1の単位が cm であるときに，次元2の単位が kg であるようなときのことをいう．ここで T を指定すると相関行列から，F を指定すると分散共分散行列から主成分分析が行われる．基本的に単位が合っていない解析の方が多いだろうから，分散共分散行列から主成分分析を行うのは限られた場合のみである．
$ rpca=prcomp(x=mtcars,scale=T)
以上のコマンドにて主成分分析が終了し，rpca には主成分スコア $x，標準偏差 $sdev (固有値の正の平方根)，$rotation (各主成分軸の固有ベクトル)，$center (解析に用いた各項目の値の平均) 等が格納される．多次元データを低次元に縮約したいだけなら，これらの項目の中では主成分スコア，$x だけを考えれば良い．以下のように打つことで，それぞれ PC1軸および PC2軸の値を確認できる．
$ rpca$[,1] $ rpca$[,2]
主成分得点のプロット
これらを以下のコマンドで散布図上にプロットする．ここで，ラベルをプロットエリア上に描くためのコマンド pointLabel() を使用するための maptools は外部パッケージなので，install.packages() で事前にインストールする必要がある．
$ library(maptools) $ par(las=1) $ plot(x=NULL,type="n",xlab="PC1",ylab="PC2",xlim=c(-5,5),ylim=c(-5,5),xaxs="i",yaxs="i",xaxt="n",yaxt="n",bty="n") $ axis(side=1,at=seq(-5,5,1),tck=1.0,lty="dotted",lwd=0.5,col="#dddddd",labels=expression(-5,-4,-3,-2,-1,0,1,2,3,4,5)) $ axis(side=2,at=seq(-5,5,1),tck=1.0,lty="dotted",lwd=0.5,col="#dddddd",labels=expression(-5,-4,-3,-2,-1,0,1,2,3,4,5)) $ points(x=rpca$x[,1],y=rpca$x[,2],pch=16,col="#ff8c00") $ pointLabel(x=rpca$x[,1],y=rpca$x[,2],labels=rownames(mtcars),cex=0.8) $ box(bty="l")
これを実行した結果は以下のようになる．Mazda RX4 と Mazda RX4 Wag のように名前が似ていて，車も似ていると思われる車が近いし，Honda Civic や Toyota Corolla のような類似する大衆車が近い位置にプロットされているのがわかる．
因子負荷量のプロット
主成分分析において新たに生じた軸，PC1 と PC2 であるが，これが何を意味するのかの解釈は観測者に委ねられる．解釈のヒントは元の軸を散布図上にプロットすることで得られる．以下のように打つ．
$ library(maptools)
$ par(las=1)
$ plot(x=NULL,type="n",xlab="PC1",ylab="PC2",xlim=c(-0.5,0.5),ylim=c(-0.5,0.5),xaxs="i",yaxs="i",xaxt="n",yaxt="n",bty="n")
$ axis(side=1,at=seq(-0.5,0.5,0.1),tck=1.0,lty="dotted",lwd=0.5,col="#dddddd",labels=expression(-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5))
$ axis(side=2,at=seq(-0.5,0.5,0.1),tck=1.0,lty="dotted",lwd=0.5,col="#dddddd",labels=expression(-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5))
$ for(i in 1:11)
$ {
$ arrows(0,0,rpca$rotation[i,1],rpca$rotation[i,2],col="#ff8c00",length=0.1)
$ }
$ pointLabel(x=rpca$rotation[,1],y=rpca$rotation[,2],labels=rownames(rpca$rotation),cex=1)
$ box(bty="l")
これを実行した結果は以下のようになる．これを確認したところ，PC1軸方向には馬力や燃費重さ等が分布している．一方で，PC2軸方向には速さとかギア数が分布している．なので，PC1軸は車のパワー，PC2軸は操作性を示す軸というように解釈できるかもしれない．
上の散布図と主成分得点の散布図を一度に簡単に描くコマンドが用意されている．以下のように打つと，同じ平面上にそれらふたつの情報をプロットできる．ただし，図の装飾ができないので以下は確認用に留め，見せる用には上のコマンドを用いてちゃんと描画すべき．
$ biplot(x=rpca)
軸の寄与率
主成分分析を用いるとこのように，高次元のデータを低次元に落とし込むことができる．このとき，上のようなたった2次元で元の11次元のどれほどの割合を説明できるのかどうかを説明する値がある．以下のようにコマンドを打つとそれを確認できる．この値を寄与率という．
$ summary(rpca)$importance
これを実行した結果は以下のようになる．2行目の Proportion of Variance が各軸の寄与率を表しており，PC1軸は，元の11次元の0.600760，すなわち約60%を説明できることとなる．3行目には累積寄与率が示されており，このとき，PC1軸とPC2軸を合わせると，全体の約84%を説明できているということになる．
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 Standard deviation 2.570681 1.628026 0.7919579 0.5192277 0.4727061 0.4599958 0.3677798 0.350573 0.2775728 0.2281128 0.1484736 Proportion of Variance 0.600760 0.240950 0.0570200 0.0245100 0.0203100 0.0192400 0.0123000 0.011170 0.0070000 0.0047300 0.0020000 Cumulative Proportion 0.600760 0.841720 0.8987300 0.9232400 0.9435600 0.9627900 0.9750900 0.986260 0.9932700 0.9980000 1.0000000
図示すると以下のようになる．
主成分得点からデータの逆算
主成分分析では，データを主成分空間上にプロットするが，これとは逆に主成分空間中のあるデータを示す点 (または任意の点) から本来のデータそのものを復元 (または新規導出) することもできる．主成分分析とは，1行a列からなるデータ X の主成分スコア S を (正規化した) X と軸の固有ベクトルA (a行1列) を掛け合わせることで求める．なので，元のデータ X はその主成分スコアにAの逆行列を掛けることで復元することができる．実際にはこの掛け算の結果を，正規化した分だけ補正する．計算は，主成分分析に使用した各変数の標準偏差 (rpca$scale) を各々の要素に掛け (対角行列化した標準偏差を掛けることと同義)，さらに各変数の平均値 (rpca$center) を加える．式で表すと以下のようになる．ここで t(A) はAの転置行列を表す．相関行列から主成分分析を行った場合，すなわち，prcomp() において scale=T で計算した場合は rpca$scale を対角行列化したものを使用する必要がある．
X= S * t(A) %*% diag(rpca$scale) + rpca$center
分散共分散行列から行った場合 (scale=F) は以下の式を考える．
X= S * t(A) + rpca$center
上で行なった PCA の結果から Mazda RX4 の値を復元する．最初に以下のコマンドで固有ベクトルAの準備をする．固有ベクトルは主成分分析の結果を格納したオブジェクト rpca の $rotation に格納されているので，それを取り出すには以下のように打つ．
$ eigen01=as.matrix(rpca$rotation[,1]) $ eigen02=as.matrix(rpca$rotation[,2]) $ eigen03=as.matrix(rpca$rotation[,3]) $ eigen04=as.matrix(rpca$rotation[,4]) $ eigen05=as.matrix(rpca$rotation[,5]) $ eigen06=as.matrix(rpca$rotation[,6]) $ eigen07=as.matrix(rpca$rotation[,7]) $ eigen08=as.matrix(rpca$rotation[,8]) $ eigen09=as.matrix(rpca$rotation[,9]) $ eigen10=as.matrix(rpca$rotation[,10]) $ eigen11=as.matrix(rpca$rotation[,11])
この情報を用いて，元のデータの1行目のデータ Mazda Rx4 の値を主成分スコアから逆算するには以下のように打つ．
$ (rpca$x[1,1]*t(eigen01) + rpca$x[1,2]*t(eigen02) + rpca$x[1,3]*t(eigen03) + rpca$x[1,4]*t(eigen04) + rpca$x[1,5]*t(eigen05) + rpca$x[1,6]*t(eigen06) + rpca$x[1,7]*t(eigen07) + rpca$x[1,8]*t(eigen08) + rpca$x[1,9]*t(eigen09) + rpca$x[1,10]*t(eigen10) + rpca$x[1,11]*t(eigen11) ) %*% diag(rpca$scale) + rpca$center
これを実行した結果は以下のようになる．元のデータを完全に再現できていることが分かる．
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 21 6 160 110 3.9 2.62 16.46 3.330669e-16 1 4 4
以上では，全11個の固有ベクトルを用いて逆算したので，完全に元のデータと一致する結果が得られたが，以下のようにすると，主成分1軸および2軸のみからデータを逆算することができる．
$ (rpca$x[1,1]*t(eigen01) + rpca$x[1,2]*t(eigen02)) %*% diag(rpca$scale) + rpca$center
これを実行した結果は以下のようになる．主成分1軸と2軸の累積寄与率が80%を超えており，高い部類なので，再現率は悪くない．
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 21.66999 5.888978 190.762 161.1851 3.949175 2.759129 16.6661 0.3380028 0.8480916 4.368929 3.729912
以上では，現実の主成分スコアからデータを逆算したが，全く架空の主成分スコアから元のデータを逆算することもできる．以下の図，赤色の点で示される主成分スコアが PC1=4，PC2=-4 からなる架空のデータを考える．
これに対して，元のデータを計算するには以下のように打つ．
$ (4*t(eigen01) + (-4)*t(eigen02)) %*% diag(rpca$scale) + rpca$center
これを実行した結果は以下のようになる．このような操作をすることで，素性のよくわかった平面上からデータを逆算できるので，売れ行きデータ等と併せることで新たな車のスペックを設計する際等に活用できるのかもしれない．
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 10.96209 8.546148 437.7041 168.9767 2.379964 5.131667 19.72803 0.2865687 -0.9197895 1.712355 1.523225