[連載]フリーソフトによるデータ解析・マイニング 第24回

Rと主成分分析

1.主成分分析とは

  観測、実験、調査では、通常個体の属性を複数の項目(変数)に分けて記録する。変数が少ない場合は、簡単なグラフや基本統計量などでデータの構造を明らかにすることができるが、変数が多くなるとデータの構造が複雑になり、解析が難しくなる。一方、変数が多くなると変数の間には相関がある可能性も増える。
  主成分分析 (principal component analysis) は、多くの変数により記述された量的データの変数間の相関を排除し、できるだけ少ない情報の損失で、少数個の無相関な合成変数に縮約して、分析を行う手法である。主成分分析の手法はホテリング (Hotelling) によって1933年頃提案された。
  変数が1つ、2つの場合は、棒グラフや散布図でデータの構造を読み取ることが可能であり、主成分分析を行う必要がないが、主成分分析の考え方を説明するため、ここでは2変数の場合の例を用いることにする。
  たとえば、表1の左側に示す3つの個体の2次元データ(x1x2 )があるとする。そのx1 を横軸、x2 を縦軸にした座標系上の散布図を図1(a)に示す。図1(a)の個体は、座標軸x1 の値x1i と座標軸 x2 の値x2ii = 1,2,3 )、つまり2次元で表記しなければならない。
  しかし、この座標系を図1(b) のようなz1z2 の座標系に変換すると1変数のみで表すことができる。

表1 二次元のデータの例
x1 x2 変換 z1 z2
個体1 1 2 2.236 0.000
個体2 2 4 4.472 0.000
個体3 3 6 6.708 0.000
分散 1 4 5.000 0.000
相関係数 1 0

図1 座標の変換

  図1 (b) の中の新しい座標z1x1x2 との関係は次の式 (合成変数、あるいは線形結合式と呼ぶ) で表すことができる。

z1 = 0.447 x1 + 0.894 x2

  この合成変数は、上記の3つの個体に限っては情報の損失なしで、表1の2次元 (x1x2 ) データを1次元 (z1 ) に縮約することができる。これは、x1x2 のピアソン相関関係が非常に強く、その相関係数が1であるからである。
  上記の合成変数z1 に個体1、個体2、個体3のx1x2 の値を代入すると、表1の右側に示すz1 の値が得られる。図1(b) に示すように、直線上に全ての点が乗ると、z2 の値は全てゼロになるので合成変数z2 は何ら情報を持っておらず、分析には役立たない。ここで言う情報とはデータのバラツキ(分散)を意味し、分散が大きいほど、情報を多く持っていると考える。
  表1で分かるように、z1 の分散はx1x2 の分散より大きい。点が完全に直線に乗る場合は情報の損失がないが、そうではない場合は情報の損失が生じる。
  その際には、より多くの情報を用いて分析を行うためには、他の合成変数(たとえばz2 )を付け加えて分析すべきである。
  上記の具体的な式を一般化するとz1= a1x1 + a2x2 となる。式の中のx1x2 は用いるデータであり、係数a1a2 は用いた既知のデータから求める。
  主成分分析では、合成変数 (線形結合式) の係数を主成分と呼ぶ。主成分分析での主な問題は、いかにデータを縮約する係数 (主成分) を求めるかである。
  通常用いる主成分分析は、データの分散が最大になるように線形結合式の係数を求める方法と、相関が最大になるように線形結合の係数を求める方法がある。
  すでに説明した2変数の線形結合式を一般化するため、m 個の個体、n 個の変数 (x1x2、・・・・・・、xn ) により構成された表2のようなデータセットがあり、これを Xm×n で表す。

表2 データ( Xm×n と記する)
x1 x2 xj xn
個体1 x11 x12 x1j  … x1n
個体2 x21 x22 x2j x2n
個体i xi1 xi2 xij xin
個体m xm1 xm2 xmj xmn

  このn 次元のデータをより低いk () 次元に縮約する線形結合の一般式を次に示す。


  この線形結合式z1 に用いる係数を第1主成分、z2 に用いた係数を第2主成分と呼ぶ。説明の便利のため、係数データを表3のように並べ、A n×k これをで示す。

表3 係数データ(A n×k と記する)
主成分
第1 第1 j  k 
x1 の係数 a11 a12 a1j a1k
x2 の係数 a21 a22 a2j a2k
xi の係数 ai1 ai2 aij aik
xn の係数 an1 an2 anj ank

  データ Xm×n と係数 An×k を線形結合式で求めた値zj を主成分得点と呼ぶ。主成分得点のデータ (Zm×k) を表4に示す。
  主成分得点 Zm×k とデータXm×n、係数 An×k との関係は、Zm×k = Xm×n An×k (行列の演算) となる。

表4 主成分得点 (Zm×k と記する)
主成分得点(Z)
z1 z2 zj zk
第1 第1 j k
個体1 z11 z12 z1j z1k
個体2 z21 z22 z2j z2k
個体i zi1 zi2 zij zik
個体m zm1 zm2 zmj zmk

  分散 (相関) を最大にする方法で主成分を求めることは、データの分散共分散行列(相関係数行列)の固有値と固有ベクトルを求める問題に帰する。その固有ベクトルが主成分であり、固有値の大小がそれに対応する固有ベクトル (主成分) に含まれる情報の多少を決める。
  データ Xm×n の分散共分散行列 (あるいは相関係数行列) の固有値は通常 λi で表記する。固有値を求めるアルゴリズムは、固有値を大小の降順に並べて返す。

  固有値は主成分得点の標準偏差の2乗に等しい。この値が大きい主成分 (固有ベクトル) ほど、元のデータの情報を多く含んでいる。最も大きい固有値に対応する主成分を第1主成分、その次に大きい固有値に対応する主成分を第2主成分と呼ぶ。
  ある主成分にデータ全体の情報がどれぐらい含まれているかは、その主成分に対応する固有値 (標準偏差の二乗) が固有値全体 (標準偏差の二乗の合計) の中でどれぐらいの割合を占めているかで説明できる。
  各固有値が固有値全体に占める割合を寄与率、その寄与率を累積したものを累積寄与率と呼ぶ。
  主成分分析を行う際には、寄与率が大きい少数個の主成分を用いる。つまり、個の主成分にデータ全体の何割の情報が含まれているかは、第主成分までの累積寄与率を用いて説明する。

2.Rの主成分分析の関数

  Rには主成分分析を行う関数 prcomp と princomp がある。ここでは、prcomp についてデータを用いて説明する。princomp の基本的な使用方法と機能は prcomp とほぼ同じである。
  主成分分析の手法を用いて、データを縮約した場合、データの構造の再現性について考察するため、人工データを用いることにする。
  たとえば、図2のように2次元平面状に、2つの同心円周上に点が散布されているとする。内側の点1〜16は半径が5cmである円周上に等間隔に点A、B、C、D、Eは半径が10cmである円弧上に等間隔に並んでいる。
  A、B、C、D、Eを観測点とし、各点1〜16までの直線距離を測ると表5のようなデータが得られる (丘本1992)。このデータを丘本の円周データと呼ぶことにする。

図2 データ cars の回帰木の折れ線


表5 丘本の円周データ
個体 A B C D E
1 50 57 74 94 112
2 57 50 57 74 94
3 74 57 50 57 74
4 94 74 57 50 57
5 112 94 74 57 50
6 128 112 94 74 57
7 140 128 112 94 74
8 147 140 128 112 94
9 150 147 140 128 112
10 147 150 147 140 128
11 140 147 150 147 140
12 128 140 147 150 147
13 112 128 140 147 150
14 94 112 128 140 147
15 74 94 112 128 140
16 57 74 94 112 128

  この例は、人工データであるためやや面白味に欠けるが、主成分分析による縮約されたデータから、元のデータ構造の再現性を説明するには都合が良い。
  ここでは、5つの変数により記述されている円周上の16個の点、円弧上の5個の点に関するデータを主成分分析で2変数に縮約した場合、円周、円弧上の点の相対位置がどの程度まで再現できるかを考察することを主な目的とする。まず次のようにデータオブジェクトを作成する。

> temp<-c( 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128, 112, 94, 74, 57, 57, 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128, 112, 94, 74, 74, 57, 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128, 112, 94, 94, 74, 57, 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128, 112, 112, 94, 74, 57, 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128)
> okamoto<-matrix(temp, 16,5, byrow=F)
> colnames(okamoto)<-c("A", "B", "C", "D", "E")

  主成分分析の関数 prcomp の使用方法は簡単であり、用いる引数も少ない。デフォルトでは、分散共分散を用いる主成分になっている。作成したデータ okamoto の主成分の要約を次に示す。

> oka.pc<-prcomp(okamoto)
> summary(oka.pc)

Importance of components:
PC1 PC2 PC3 PC4 PC5
Standard deviation
68.0019 40.8155 3.83644 1.30318 0.30578
Proportion of Variance
0.7332 0.2641 0.00233 0.00027 0.00001
Cumulative Proportion
0.7332 0.9974 0.99972 0.99999 1.00000

  ソフトによっては、主成分分析の結果として固有値を返すが prcomp ではその代わりに標準偏差 (standard deviation) を返す。
  標準偏差は、各主成分得点の標準偏差で、固有値の正の平方根に等しい。元の変数がn 個あると、非ゼロである固有値および主成分(固有ベクトル)はk () 個求まる。
  主成分分析の目的は、できるだけ少ない主成分に、もとの変数の情報を吸収することにある。
  主成分分析を行う際には、用いる第q 番目の主成分までに、もとの変数の情報がどれだけ吸収できたかが問題となる。その判断のために、寄与率と累積寄与率に関する情報が必要となる。
  返された結果の中の Proportion of Variance は、各標準偏差の二乗が標準偏差の二乗の合計に占める割合(寄与率)であり、Cumulative Proportion は累積寄与率である。
  返された結果から分かるように第2主成分までの累積寄与率が0.997(99,7%)であり、5次元データの情報のほとんどが第1、2主成分に縮約されている。
  関数 prcomp で計算された主成分(固有ベクトル)は $rotation に、列を単位に記録する。主成分は、変数の相互関係を分析するのに用いる。丘本円周データの主成分を次に示す。

> oka.pc$rotation

PC1 PC2 PC3 PC4 PC5
A 0.3648367 6.201120e-01 0.5837693 -3.397958e-01 -0.1615786
B 0.4805596 3.397958e-01 -0.1641399 6.201120e-01 0.4920575
C 0.5214531 -8.326673e-17 -0.5143375 1.110223e-15 -0.6808404
D 0.4805596 -3.397958e-01 -0.1641399 -6.201120e-01 0.4920575
E 0.3648367 -6.201120e-01 0.5837693 3.397958e-01 -0.1615786

  主成分分析では、情報がより多く含まれている上位のいくつかの主成分を用いる。分析にはそれらの棒グラフや散布図が多く用いられている。
  円周データでは第2主成分までの累積寄与率が99.7%に上るので、図3に第1、2主成分の散布図のみを示す。

> plot(oka.pc$rotation[,1],oka.pc$rotation[,2],type="n")
> text(oka.pc$rotation[,1],oka.pc$rotation[,2],colnames(okamoto))

  この図3を、図2と比べると、主成分の散布図円弧が若干変形されたものの、図3を、横軸を軸として180度回転するとA、B、C、D、Eの相対位置は元の点の位置とほぼ一致し、元のデータの相対的な構造が再現されていると言える。

図3 第1,2主成分の散布図

  主成分が、線形結合式の合成変数を求める係数であるので、用いたデータと主成分の結果で主成分得点を求めることは可能であるが、通常の主成分分析のソフトパッケージは、主成分得点を返す機能を揃えている。関数 prcomp では、主成分の得点は、$x に列単位で記録する。データ okamoto の主成分得点を次に示す。主成分得点は、個体のデータ構造に関する情報を縮約したものである。

> oka.pc$x

PC1 PC2 PC3 PC4 PC5
[1,]
-65.348592 -5.101938e+01 -2.879591 -1.876805e+00 -0.13974065
[2,]
-91.201607 -3.109924e+01 3.874461 -2.310244e+00 -0.07364076
[3,]
-100.751885 3.197442e-14 7.364914 1.598721e-14 0.25640318
[4,]
-91.201607 3.109924e+01 3.874461 2.310244e+00 -0.07364076
[5,]
-65.348592 5.101938e+01 -2.879591 1.876805e+00 -0.13974065
[6,]
-29.708700 5.694019e+01 -5.484544 -5.612449e-01 -0.25084539
[7,]
7.557865 5.248045e+01 -3.722346 -1.342714e+00 0.52231640
[8,]
40.168493 4.238022e+01 -1.114172 -6.460407e-01 0.02797145
[9,]
65.140371 3.002037e+01 1.197715 -1.130112e+00 -0.21794250
[10,]
80.741814 1.518008e+01 2.724255 -2.550000e-01 0.29651471
[11,]
86.052595 -3.552714e-14 3.443529 -7.105427e-15 -0.58566967
[12,]
80.741814 -1.518008e+01 2.724255 2.550000e-01 0.29651471
[13,]
65.140371 -3.002037e+01 1.197715 1.130112e+00 -0.21794250
[14,]
40.168493 -4.238022e+01 -1.114172 6.460407e-01 0.02797145
[15,]
7.557865 -5.248045e+01 -3.722346 1.342714e+00 0.52231640
[16,]
-29.708700 -5.694019e+01 -5.484544 5.612449e-01 -0.25084539

  さて、各個体(図2に示された16個の点)の構造を合成変数で縮約した主成分得点で再現するため、次のように第2主成分までの主成分得点の散布図を作成してみる。

> plot(oka.pc$x[,1],oka.pc$x[,2],type="n")
> text(oka.pc$x[,1],oka.pc$x[,2],1:16)

図4 第1,2主成分得点散布図

  図4で分かるように、縦軸を軸とし180度回転すると円形が若干変形されたものの、16個の点の相対位置は図2と変わらない。
  ここでは5変数 (次元) のデータを2変数に縮約している。縮約した2次元のデータは、多少歪みはあるものの、元のデータ構造をある程度再現している。
  主成分分析は、変数が多いとき情報の損失を最小限に押さえながら少ない合成変数に縮約する方法であるため、元のデータ構造が100%正確に再現できないことと、再現されているのは元の変数、または個体間の相対的な関係に過ぎないことを強調しておきたい。
  主成分分析を行うとき、いくつの合成変数 (主成分) を用いて分析するのが適切であるかに関しては明確な決まりがない。
  分散共分散行列を用いる場合は、一般的には累積寄与率70%〜80%を大まかな目安とし、累積寄与率がこれを超える主成分まで用いて分析をすることが多い。
  相関行列を用いた主成分分析の場合は、固有値の値が1前後になる主成分まで用いるのが1つの目安である。
  関数 prcomp には引数 scale があり、データの標準化 (データのスケールを統一) が必要なときは scale = TRUE を指定する。デフォルトには scale = FALES になっている。scale = TRUE にすると元のデータの相関行列を用いた主成分に等しい。これは関数 princomp の引数 cor = TRUE の結果と対応する。
  関数 prcomp を用いた主成分分析では、関数 predict を用いて、あるデータに基づいて作成した合成変数に、新しいデータを当てることが出来る。その書き式を次に示す。

predict(object,newdata,..)

  たとえば、丘本の円周データの個体1〜10を用いて、合成変数を作成し、残りの11〜16の個体を、求めた合成関数に当てた2次元の主成分得点の散布図を図5に示す。

> oka.pc1<-prcomp(okamoto[1:10,])
> oka.pc2<-predict(oka.pc1,okamoto[11:16,])
> plot(oka.pc1$x[,1],oka.pc1$x[,2],type="n")
> text(oka.pc1$x[,1],oka.pc1$x[,2],1:10)
> plot(oka.pc2[,1],oka.pc2[,2],type="n")
> text(oka.pc2[,1],oka.pc2[,2],11:16)

図5 predict の結果の散布図

  関数 prcomp の結果については、主成分と主成分得点を同一の画面上で散布図を作成する biplot 関数を用いることができる。関数 biplot を用いた結果を次に示す。

> biplot(oka.pc)

図6 関数 biplot による散布図

  主成分分析で求まる主成分および主成分得点の正負の符号は、固有値及び固有ベクトルを求めるときのアルゴリズムが異なると逆になる場合がある。
  つまり、個体の散布図を描いたときに、異なるアルゴリズムによる結果の上下が逆になったり、左右が逆になったりすることがある。主成分分析で行う分析は、変数間、個体間の絶対的関係ではなく、相対的関係であるため、分析には問題がない。

引用文献:
丘本正(1991): 多変量解析の諸方法のモデル再現性に関する数値実験、行動計量学、Vol.18、No.2、P.47-56.