Ｒと主成分分析
|
１. 主成分分析とは |
観測、実験、調査では、通常個体の属性を複数の項目（変数）に分けて記録する。変数が少ない場合は、簡単なグラフや基本統計量などでデータの構造を明らかにすることができるが、変数が多くなるとデータの構造が複雑になり、解析が難しくなる。一方、変数が多くなると変数の間には相関がある可能性も増える。
主成分分析(principal component analysis)は、多くの変数により記述された量的データの変数間の相関を排除し、できるだけ少ない情報の損失で、少数個の無相関な合成変数に縮約して、分析を行う手法である。主成分分析の手法はホテリング(Hotelling)によって1933年頃提案された。
変数が１つ、2つの場合は、棒グラフや散布図でデータの構造を読み取ることが可能であり、主成分分析を行う必要がないが、主成分分析の考え方を説明するため、ここでは2変数の場合の例を用いることにする。
たとえば、表1の左側に示す3つの個体の2次元データ(
、
)があるとする。そのを横軸、を縦軸にした座標系上の散布図を図１(a)に示す。図１(a)の個体は、座標軸の値
と座標軸の値
（
)、つまり2次元で表記しなければならない。
しかし、この座標系を図1(b)のような
、
の座標系に変換すると1変数のみで表すことができる。
|
表1 二次元のデータの例 |
|||||
|
|
|
|
変換 |
|
|
|
個体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 座標の変換 |
||
|
|
変換 |
|
図１(b)の中の新しい座標 と
、 との関係は次の式(合成変数、あるいは線形結合式と呼ぶ)で表すことができる。
=0.447+0.894
この合成変数は、上記の3つの個体に限っては情報の損失なしで、表1の2次元(
、 )データを1次元( )に縮約することができる。これは、 と
のピアソン相関関係が非常に強く、その相関係数が1であるからである。
上記の合成変数 に個体１、個体２、個体３の 、
の値を代入すると、表１の右側に示す の値が得られる。図１(b)に示すように、直線上に全ての点が乗ると、 の値は全てゼロになるので合成変数は何ら情報を持っておらず、分析には役立たない。ここで言う情報とはデータのバラツキ（分散）を意味し、分散が大きいほど、情報を多く持っていると考える。
表１で分かるように、 の分散は 、
の分散より大きい。点が完全に直線に乗る場合は情報の損失がないが、そうではない場合は情報の損失が生じる。
その際には、より多くの情報を用いて分析を行うためには、他の合成変数(たとえば )を付け加えて分析すべきである。
上記の具体的な式を一般化すると =
+ 
となる。式の中の
、
は用いるデータであり、係数
、
は用いた既知のデータから求める。
主成分分析では、合成変数(線形結合式)の係数を主成分と呼ぶ。主成分分析での主な問題は、いかにデータを縮約する係数(主成分)を求めるかである。
通常用いる主成分分析は、データの分散が最大になるように線形結合式の係数を求める方法と、相関が最大になるように線形結合の係数を求める方法がある。
すでに説明した2変数の線形結合式を一般化するため、
個の個体、 
個の変数( 、 、･･････、
)により構成された表2のようなデータセットがあり、これを
で表す。
|
表2 データ( |
|
|
この次元のデータをより低い
(
)次元に縮約する線形結合の一般式を次に示す。
この線形結合式 に用いる係数を第1主成分、 に用いた係数を第2主成分と呼ぶ。説明の便利のため、係数データを表3のように並べ、これを
で示す。
|
表3 係数データ( |
|
|
データと係数を線形結合式で求めた値
を主成分得点と呼ぶ。主成分得点のデータ(
)を表4に示す。
主成分得点とデータ、係数 との関係は、
= (行列の演算)となる。
|
表4 主成分得点( |
|
|
分散(相関)を最大にする方法で主成分を求めることは、データの分散共分散行列(相関係数行列)の固有値と固有ベクトルを求める問題に帰する。その固有ベクトルが主成分であり、固有値の大小がそれに対応する固有ベクトル(主成分)に含まれる情報の多少を決める。
データの分散共分散行列(あるいは相関係数行列)の固有値は通常
で表記する。固有値を求めるアルゴリズムは、固有値を大小の降順に並べて返す。
固有値は主成分得点の標準偏差の2乗に等しい。この値が大きい主成分(固有ベクトル)ほど、元のデータの情報を多く含んでいる。最も大きい固有値に対応する主成分を第１主成分、その次に大きい固有値に対応する主成分を第２主成分と呼ぶ。
ある主成分にデータ全体の情報がどれぐらい含まれているかは、その主成分に対応する固有値(標準偏差の二乗)が固有値全体(標準偏差の二乗の合計)の中でどれぐらいの割合を占めているかで説明できる。
各固有値が固有値全体に占める割合を寄与率、その寄与率を累積したものを累積寄与率と呼ぶ。
主成分分析を行う際には、寄与率が大きい少数個の主成分を用いる。つまり、
個の主成分にデータ全体の何割の情報が含まれているかは、第主成分までの累積寄与率を用いて説明する。
2．Rの主成分分析の関数
|
Rには主成分分析を行う関数prcompとprincompがある。ここでは、prcompについてデータを用いて説明する。princompの基本的な使用方法と機能はprcompとほぼ同じである。
主成分分析の手法を用いて、データを縮約した場合、データの構造の再現性について考察するため、人工データを用いることにする。
たとえば、図２のように２次元平面状に、２つの同心円周上に点が散布されているとする。内側の点1〜16は半径が5cmである円周上に等間隔に点A、B、C、D、Eは半径が10cmである円弧上に等間隔に並んでいる。
A、B、C、D、Eを観測点とし、各点1〜16までの直線距離を測ると表5のようなデータが得られる（丘本1992）。このデータを丘本の円周データと呼ぶことにする。
|
図2 丘本の円周データ図 |
|
|
|
表5 丘本の円周データ |
|
|
この例は、人工データであるためやや面白味に欠けるが、主成分分析による縮約されたデータから、元のデータ構造の再現性を説明するには都合が良い。
ここでは、５つの変数により記述されている円周上の16個の点、円弧上の5個の点に関するデータを主成分分析で2変数に縮約した場合、円周、円弧上の点の相対位置がどの程度まで再現できるかを考察することを主な目的とする。まず次のようにデータオブジェクトを作成する。
|
> temp<-c( |
|
> okamoto<-matrix(temp, 16,5, byrow=F) |
|
> colnames(okamoto)<-c("A", "B", "C", "D", "E") |
主成分分析の関数prcompの使用方法は簡単であり、用いる引数も少ない。デフォルトでは、分散共分散を用いる主成分になっている。作成したデータokamotoの主成分の要約を次に示す。
|
>oka.pc<-prcomp(okamoto) |
|
>summary(oka.pc) |
ソフトによっては、主成分分析の結果として固有値を返すがprcompではその代わりに標準偏差(standard
deviation)を返す。
標準偏差は、各主成分得点の標準偏差で、固有値の正の平方根に等しい。元の変数がn個あると、非ゼロである固有値および主成分(固有ベクトル)は
( ) 個求まる。
主成分分析の目的は、できるだけ少ない主成分に、もとの変数の情報を吸収することにある。
主成分分析を行う際には、用いる第q番目の主成分までに、もとの変数の情報がどれだけ吸収できたかが問題となる。その判断のために、寄与率と累積寄与率に関する情報が必要となる。
返された結果の中のProportion of Varianceは、各標準偏差の二乗が標準偏差の二乗の合計に占める割合(寄与率)であり、Cumulative
Proportionは累積寄与率である。
返された結果から分かるように第2主成分までの累積寄与率が0.997(99,7%)であり、5次元データの情報のほとんどが第1、2主成分に縮約されている。
関数prcompで計算された主成分(固有ベクトル)は$rotationに、列を単位に記録する。主成分は、変数の相互関係を分析するのに用いる。丘本円周データの主成分を次に示す。
|
>oka.pc$rotation |
主成分分析では、情報がより多く含まれている上位のいくつかの主成分を用いる。分析にはそれらの棒グラフや散布図が多く用いられている。
円周データでは第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 |
|
|
さて、各個体（図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%を大まかな目安とし、累積寄与率がこれを超える主成分まで用いて分析をすることが多い。
相関行列を用いた主成分分析の場合は、固有値の値が１前後になる主成分まで用いるのが１つの目安である。
関数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.