 |
|
| Chemometrix > ケモメトリックス > PCA | |
 |
7. 主成分分析(Principal Component Analysis) |
||
| 主成分分析(PCA)では、２つの重要な仕事を行えます。最初の１つは、解析結果を視覚化して探索的データ解析を行うことです。PCAはHCAのデンドログラムように、サンプル間または変数間の関係をグラフとして表示できます。2番目は、データの次元数を減らす方法を提供します。PCAは、元の説明変数の変数の最大の変動量を占めるの直線の組み合わせを見つけます。多変量距離がHCAに対して重要なのと同様に、分散の概念は、主成分分析に対して重要です。 n個の要素を含むベクトルxの場合、ベクトルに対応する変動(variation)は一般的に以下のように定義されます  この定義は、下記のように定義される数学的な分散(variance)と対比させられます｡  もし、データを平均化またはオートスケールで常に前処理しているのであれば、変動量と分散の区別について議論する余地があります。これら両方の前処理オプションはそれぞれの変数の平均を引きますので、2つの質はn-1の係数だけ異なっています。しかし、データを前処理しない場合、またはレンジスケールか分散スケールを選択している場合、変動量と分散の違いに気が付くでしょう。今後の説明で、分散という用語は式２で計算されている量を表現するのに使用します。 前処理の説明については、「前処理」を参照してください。 HCAが、多変量距離の類似性にもとづいて計算されるのと同様に、PCAでは変動量の情報にもとづいて計算します。変動量は、適切または不適切として分類されます。ノイズは不適切な変動量の例です。 大きなデータセットを視覚化することが便利という点は納得できると思いますが、なぜ次元を減らすことが有効なのか不思議に思われるでしょう。例えば、複数のセンサーによって測定された多変量なデータセットを考えます。この複数のセンサー(変数)が果物の新鮮さに応答します。しかし、複数のセンサーで得られた情報には、部分的に重複する場合があります。統計的な表現をすると、それらの測定された変数には相関があります。複数のサンプルのデータが、50個の部分的に重複するセンサーから収集される場合、実測の分散パターンはおそらく50の独立な現象によるのではなく、むしろそれより小さい数のファクターによって要約できます。この数は、データの固有の次元として表現されます。PCAでは、これらのファクターを見つけ出す方法を提供します。 PCAによって発見される元の説明変数の直線的な組合わせには、ファクターや因子、主成分、ローディング、固有ベクトル等の用語が使用されます。PCAの固有ベクトルには、いくつかの望ましい特性があります。第一に、互いに直交する性質があります｡ つまり、相関がありません｡ 第二に、固有ベクトルは分散を減少する方向で計算される点です。よって、最初の固有ベクトルが分散の最大値になり、その後の固有ベクトルは、残りのデータの分散を少なくするように計算されます。固有ベクトルの数は、元の変数やサンプルの少ない方の数になります。 相互関係のない変数の直線の組合わせ(即ち固有ベクトル)を見つけることで、変数間の相関性によって生じる困難を避けられます。 PCAの結果を視覚化する前に、重要な変数を座標軸に使用した2次元または3次元の散布図をいくつか作成します。この方法により、サンプル間の重要な関係を識別できます。しかしながら、ここで重複があると、どの変数が最も重要なのか前もって知ることができません。さらに、座標軸の変数の組み合せの数は、変数が増大すると大幅に増えます。ここで、PCAの結果の最初の2つか3つの主成分を、元の変数の代わりの座標軸とすれば、表示されるサンプル間の関係は得られた分散の観点から最適化されます。これは視覚化の技術を使用しての主成分分析の長所です。 主成分分析では、多変量データセットの説明変数ブロックにおける変化を可能な限り最善の表示で提供します。この表示により、データのクラスタリングやアウトライヤーのサンプルが明確になります。また、主成分分析で明確になるデータのパターンに化学的(または生物学あるいは物理学)意味を与えるが可能です。そして、測定箇所のどの部分がノイズであるかを推定できます。 備考:PCAは、SiMCAやPCRの基礎となっているため、PCAの基本的な理解は非常に重要です。 最大の分散の軸の考えを使用して、バナナの次元数を減らすことを考えてみます。バナナの長さの軸が、最初の主成分(PC1)になります｡ 最初の座標軸に対して垂直な2番目の座標軸は、2番目の主成分(PC2)を意味します。第2主成分はその果房の先から得られる高さです。3番目の主成分は第1主成分と第2主成分の両方に対し直角に引かれた座標軸で、バナナの厚さになります(PC3)。図17は、主成分軸での表現に基づく、3つの異なる二次元表示の遠近図を示しています。図18 は3つの主成分の座標軸を重ねた遠近表示です。 図18 で示されているバナナは、単に3つの新しいプロットの座標軸に変換されただけです。この場合、次元数の縮小は行なっていません｡ x、y、z座標をもつバナナは､ 3つの次元まで表現できます｡ しかし､ もし、3次元のバナナを平面表示に縮小する必要がある場合(例えば、影として表現する場合)、 図17の左下に示されている第1主成分と第2主成分の表示が、最大限にバナナを表現していることになります。  図 17 バナナの表示  図18 3つの主成分軸で表現したバナナ |
||
1. 数学的背景 |
||
PCAは、スコア行列Tとローディング行列Lの転置行列の2つから、Xを表現する考えに基づいています。 [式 14] X行列の大きさは、n行とm列があり、それぞれサンプルの数と説明変数の数に対応しています。スコア行列Tにはn行とg列、ローディング行列Lにはm行とg列、Lの転置行列にはg行とm列があります。ここで、gの値は、説明変数の数とサンプル数の小さいほうの数です。 スコア行列とローディング行列の最初のk列が(ここでkはgより小さい)使用された場合、近似の  が計算できます。 [式 15] ここで､  はXの概算値で、データ行列の次元が縮小されたと考えます。(実際に減少したのはTとLです。)PCAの視覚化 Lの各列は主成分を表わし、元の変数の直線的な組み合せの新しいファクターです。Lの最初の列のm個の要素を含むローディングは、元のそれぞれの変数が最初の主成分でどのように影響されるかを示しています。スコア行列Tは、ローディングによって定義された座標軸へのサンプル投影です。それぞれのサンプルは、個々の新しい座標軸に座標を持ちます。Tの列にはこれらの座標が示されます。 次の図は、6つのサンプルと2つの変数を含むデータセットについて、PCAの概念を説明しています。図19 の2つのファクター(PC1とPC2)が、元の変数(var1とvar2)の座標軸上に示されています。PC1がデータの主の広がりを表現している点に注意してください。  図19 6点(白丸)の主成分軸の図解 図20は、2つのサンプルのスコア、(tij) が、元の変数の空間でのデータ位置を、主成分の観点からどのように表現されているかがわかるように点線で記されています。smp1はPC1とPC2のスコアがほぼ等しいですが、smp2ではPC2よりPC1でスコアがずっと大きくなっています。  図20 スコア: tikはk番目の主成分軸上のi番目のサンプルのスコア ローディングは、それぞれの主成分軸に関する元の変数の相対的な貢献を表わしています｡ 1つの説明変数から成っているファクターはまれで、それぞれのファクターは、元の説明変数全ての(またはほとんどの)直線的な組み合せになります。 得られるファクターに対する変数の貢献度については、そのファクター座標軸における任意のポイントを、元の変数の座標軸に投影して示されます。図21から、var2はvar1より、PC1に寄与しているのが明らかです｡ 変数のローディングは必ず1と-1の間の値になります。そしてその絶対値は、その変数がどの位ファクターに負荷しているかの尺度になります。1の値のローディングは、その変数が主成分の軸と一致していることを意味します。  図21 ローディング: lkjはk番目の主成分軸上のj番目の変数のローディング PC1は、それぞれのサンプルポイントと、その線までの距離の二乗の合計を最小にする線です。したがって、最初のファクターの方向はデータの一番大きい広がりの線になります。2番目のファクターの座標軸(PC2)は、エルミット行列の定義から最初のファクター(PC1)に対して垂直になります。その方向はデータの残っている変動を表現する、もう1つの線になります。 それぞれのファクターの座標軸に関連するのが、そのファクターによってとらえられた分散の量を示す固有値です。固有値は、ファクター数の増加と共に減少します。表現されていない残りの全分散の合計が小さくなるほど、固有値も小さくなります。 PCAによるモデルの構築 情報を著しく損なう事なく、後の複数のファクターの省略を認めるかは､ 最初のいくつかの主成分における分散の割合に依存します｡ データ行列の説明変数のブロックは、残差エラー、E を含んで再構築されます。  [式 16] 下付きの記号は、スコアやローディング行列が、kまでのファクターを含んで縮小されているだけでなく、残差行列がkに依存していることも示しています。 同様に重要な利点は、予測の主成分分析モデルを作成することです。縮小されたローディング行列Lkは、新しいデータを解析するために後で使用できます。 ある特定の原材料を信頼できる提供者から得ていたのですが、後になってその提供者を変更するとします。以前の提供者からの原材料について、異なるバッチのスペクトルデータからの主成分分析モデルを作成し、そこから新しい材料のスペクトルを収集すると、新しい材料と古い原材料の比較ができます。それぞれの新しいサンプルのスペクトルのスコアを以下の式から予測できます。  [式 17] これは、PCAモデルで定義された空間へ新しいデータの投影として説明できます。もし、予測のスコアが、許容される既知サンプルのスコアの近くに存在していれば、スペクトルに大きな違いはありません。 モデルのバリデーション モデルの構築には、いくつかの項目について慎重に検討する必要があります。
ファクターの最適数を予測する ファクターを追加するのを、いつストップするかは、非常に重大な問題です。この決定を自動的にする多くの方法が報告されています）。「PCR/PLS」では、2種類のストップの判定基準を説明します｡ しかし、自動ストップ方式は完璧ではないため、いくつのファクターを使用するかを決める前に、必ず異なったkの設定で複数の診断的方法を慎重に検討してください。 サンプル残差 サンプル残差の分散は、直接に残差行列Ekから導かれます。表記をより容易にするために、下付きのkの記号をなくしハットの付いたシンボルがkファクターによる値を示します。  のi番目の行のベクトル( )は、サンプルの元のデータ( )とそのkファクターによる予測値( )の差です。 [式 18] ここで、サンプル残差の分散は、  [式 19] 全体の分散は、トレーニングセット全体として以下のように計算できます:  [式 20] あるサンプルのサンプル残差がs0より大きい場合、そのサンプルがアウトライヤーと考えられます。つまり、それがトレーニングセットの他のサンプルと同じ母集団に属さない場合です。この場合、F検定が、2つの分散が大きく異なるかどうかを決めるのに使用されます。  [式21] ここで、n/(n-k-1)は､ 考慮の対象になるサンプルがs0の計算に含まれているために加えられています。もし、式21の左辺が、(指定された確率の1とn-k-1の自由度に基づいて)F分布表からの得られる値と同じに設定されている場合、判定値のサンプル残差は式21を変形して求められます。  [式 22] これが、後でサンプル残差が「大きすぎる」かどうかを判定する基準になります。もし、サンプル残差がscritを越えている場合、そのサンプルはアウトライヤーの可能性があります。 マハラノビス距離 マハラノビス距離は、k個のファクターの主成分空間におけるサンプルの平均位置からの1つのサンプルの距離です。  [式23] ここで、Sはスコアの分散共分散行列です。 マハラノビス距離は普通分布していると考えられ、アウトライヤーの判定値は、自由度kのカイ二乗分布から求められます。 エラーの二乗和 式19のサンプル残差の分散は、それぞれのサンプルについての要素を持っている列ベクトルです。もし、個々の要素が二乗され、合計が自由度で修正されれば、予測残差二乗和(predicted residual error sum of squares(PRESS))を計算できます。  [式24] ファクターの数が増加するに従い、近似の質が高くなり、PRESSはゼロに近づいていきます。 バリデートされたエラーの二乗和 式19で説明されたサンプル残差の分散が、クロスバリデーションの際に増加した場合、それに対応するバリデートされた予測残差二乗和(VPRESS)が計算されます。  [式25] 式25で定義された式には、混乱を避けるために下付き記号のcvが加えられています。ファクター数が減少すると、VRESSは常に単調に減少しません。最小に達すると、その後増加する場合があります。つまり、ファクターをそれ以上加えても、モデルの予測能力は改善されません。 サンプル残差とマハラノビス距離は、サンプルに対しての計算方法です。モデリングパワーは主成分数によって変化し、変数に対応しています｡ 一般に、モデリングパワーは使用するファクターの最適数を決める助けにはなりませんが、重要な変数を指摘してくれます。この量に関して最初に、残差行列Eのj番目の列を使って、変数の残差分散、  を計算します。モデリングパワー  [式26] この変数の全分散(  )は、 [式 27] モデリングパワー(  )は、以下のように定義されます。 [式 28] データのモデル情報に対する変数の影響が増加すると、モデリングパワーは1に近づきます。モデリングパワーはしばしば負の値になります。ランダムなデータでも、「重要な」変数が高いモデリングパワーを示しますので、絶対的な判定基準は特定されません。 主成分分析での予測 主成分分析での予測は、トレーニングセットのモデルのローディングによって定義された主成分空間に未知サンプルを投影して行われます。その空間に位置している未知サンプルについての判断は、「サンプル残差」で概説されている方法に似たやり方で、残差の大きさに基づいて決められます。 はじめに、未知サンプル(つまり、トレーニングセットにないサンプル)を考える時、判断となるサンプル残差、scritの表現は以下のようになります。  [式 -29] 次に、図3-22で示されているような、AとBの2つの未知サンプルが、トレーニングセットの2つの主成分によって定義された空間に投影されている例で説明します。  図 22 未知サンプルAとBに関する2つの主成分モデル 未知サンプルAは明らかにトレーニングセットの領域に位置していますが、Bは異なります。しかし、両方ともscritより少ない、(式19からの平方根から計算された)サンプル残差と考えます｡ このような状況を処理するために、未知のサンプルの残差の計算を次項で説明されているスコアハイパーボックスに基づいて修正します。図23は、そのハイパーボックスを示しています。 図23 図22のハイパーボックス スコアのハイパーボックス それぞれの主成分に関しての、スコアのハイパーボックスの下限は、  [式30] そして、上限は  [式31] ここで、tminとtmaxはスコアのそれぞれの最小値と最大値で、stはスコアの標準偏差です｡ t*は、トレーニングセットのサンプル数に基づく自由度としてt分布から得られます。そして、cは標準偏差マルチプライヤーと呼ばれる定数です。c値の範囲は変化できるため、負または正の値を指定してtminとtmaxで示されるの境界を縮小または拡大できます。二次元のスコアプロットとして表示すると、ハイパーボックスと呼ばれる長方形の境界が形成されます。 予測におけるサンプル残差の増加 ポイントBは、明らかに図23で示されているハイパーボックスの外側に位置しています。サンプルからハイパーボックスの端までの距離に比例してサンプル残差の量が増大すると、Bがアウトライヤーとして判別される材料が増えます。増大させられたサンプル残差(距離とも呼ばれます)は以下の式から計算されます。  [式32] シグマはk主成分までの加算です。ジグマ記号の前の比で、単位を同じにします。 NiPALS (Nonlinear iterative Partial Least Squares) 主成分分析の計算には、NiPALSを使用する場合があります。なぜなら、NiPALSでは、可能な全てのファクターを計算しないで、最初のk主成分を見つけられます(実際に、全てのファクターを計算するのは非常に効率が悪いです。 |
||
2. PCAのモデル構築 |
||
|
PCAの計算が実行されると、多くの計算結果を表示できます。PCAによって計算されたオブジェクトは、サンプルのアウトライヤーを見つけたり、ファクターの適切な数を選んだり、変数を除外する等の判断をする際に役立ちます。それらをどのように調査するかの考えを含めて説明します。モデルの構築は反復作業である点を、心に留めておいてください。主成分分析の計算を一回だけ実行し、すぐに予測を行うのはまれです。つまり、代わりに時間のほとんどをモデルの適正化に使用して、サンプルと変数、環境設定パラメータの「最適」な設定を見つけることです。 固有値 固有値オブジェクトのテーブルを、表2に示します。Varianceの名前が表示されている列には、それぞれの主成分に対応する固有値が含まれています。これらは主成分によって表現された変動量です。Percentが先頭の列は、変動量をパーセント、Cumulativeの列はパーセントの累積(寄与率)を表わしています。一番大きい分散は、常に第1主成分に関連しています。より多くのファクターを使用すると、関連している分散は減少します。この場合、第1主成分が全体の分散の約53%を占めています。そして第5主成分までで95%の分散が占められます。 表2 固有値のテーブル表示  エラー PCAのエラーとは、式24で定義されているPRESS(予測残差平方和)です。ここでは、この測定を、Press Calと呼び、図24で示しています。  図24 Press Calのラインプロット クロスバリデーションをPCAの計算で行った場合、エラーには、式25で定義された量も含みます。ここではPress Valと呼びます。図25では、Press CalとPress Valがラインプロットとして重ねて表示されています。 Press ValはPress Calよりも大きい値になっています。  図25 Press CalとPress Valのラインプロット スコア スコアはサンプル間の関係を示すため、探索的データ解析に不可欠です。スコアプロットを調べる際には、必ずその調査の目的を心にしっかり留めておいてください。もし、最終目的が既知のクラスに対応するプロットの分類とサンプルのグルーピングの場合、良いKNNやSiMCAのモデルを構築できる可能性は高いです。しかし、回帰が最終目的であれば、複数のサンプルが単一の均質な形になるのが望ましいです。この場合、アウトライヤーや複数のサンプルによってクラスターが2つ以上存在すると、モデルの質を落とす可能性があります。もし十分な大きさのクラスターが複数個存在する場合、それぞれのグループにおいて個々に回帰モデルを作成します。三次元のスコアプロットを、下図に示します。  図26 PCAスコアの三次元プロット ローディング 三次元散布図のローディングプロットを、図27に示します。  図27 PCAローディングの三次元プロット ローディングは各主成分に対し、どの変数が最も重要で、どれがほとんど役に立たないかを指摘します。元のデータとローディングのラインプロットを比較すると、特定の主成分によってとらえられているデータ特長がわかります。このような比較を、図28に示します。  図28 生データのラインプロットとローディングベクトルの比較 X残差 k個までのファクターに基づく元の説明変数のブロックの概要(  )は、以下の式から計算されます。 [式 33] Xの残差  は元のデータとk個のファクターからの概要の差です。 [式 34] Xの残差は、kファクターモデルによって、比較的適合の一致しないサンプルと変数を指摘する際に有効です。また、元のデータと近似されたデータとの差が、正当な大きさであるかどうかを確かめられます。X残差をラインプロットとして表示した例を下図に示します。  図 29 X残差 もし、オートスケールや分散スケール、レンジスケールを前処理として選択した場合、X残差は元のデータと同じ単位になりません。よって、2つの大きさの比較に、意味はありません。しかし、前処理を行なわない場合や、平均化の場合、この比較は有効かつ有益です。 アウトライヤーの診断 k個のファクターによって、各サンプルがどの程度うまく近似されているかを評価する、2つの測定(サンプル残差とマハラノビス距離)方法について説明します。これらの2つの測定を軸にした二次元のプロットにした例を図29に示します。  図30 アウトライヤーの診断 スレッシュホールドの1つまたは両方の外側に位置するサンプルは、アウトライヤーの可能性のあるサンプルです。この場合、サンプル残差のスレッシュホールドは95%の確率限界に基づいていますので、「正常」サンプルの5%は、そのカットオフ値の外側になります。よって、大きなデータセットの場合、たった1つのスレッシュホールドをわずかだけ越えているサンプルは、多分「正常」なサンプルです。しかし、1つのスレッシュホールドを大幅に越えているか、または両方のカットオフ値を越えているサンプルは、アウトライヤーな可能性が大きいです。したがって、必ずこれらのサンプルを詳しく調べて、他のサンプルとどのように違うのかを理解してから、これらのサンプルを除外して主成分分析を再実行することを考えます。 |
||
３．PCAでの予測結果 |
||
|
PCAを実行すると、モデルが構築されます。PCAを実行した後で、 予測を行うには、モデルとターゲットになるデータが必要です。つまり、1つ以上のサンプルから成りたつXブロックを含むデータセットが必要です。このデータセットのXブロックには、モデルが作成されたデータセットと同じ説明変数と数を含んでいる必要があります。予測されるデータに、目的変数やクラス変数を含ませる必要はありません。 スコア 予測が実行されるとスコアが計算されます。また、スコアにはハイパーボックスを規定できます。ボックスの大きさは、式30と式31で設定されている2つのパラメータの設定に従っています。「スコアハイパーボックス」を参照してください。二次元散布図のプロットで描かれたスコアハイパーボックスを、図31に示します。  図31 ハイパーボックスを含むスコアの表示 予測されたスコアがスレッショホールドの外側に位置する場合、その新しいサンプルはトレーニングセットのサンプルと異なることを示唆します。ハイパーボックスを大きく越えているサンプルに関しては、X残差を必ず確かめます。 X残差 図29と図32を、y座標軸のスケールの違いに特に注意して比較してください。トレーニングセットと大幅に異なるサンプルは、下図のように非常に大きな残差を示します。  図32 予測におけるX残差 予測時のアウトライヤー診断 下図には、式32で定義された増加したサンプル残差と式29で定義されたスレッシュホールドを含んでいます。x軸がサンプル番号、y軸がサンプル残差の二次元散布図のプロットとして示されます。ここで、スレッシュホールドのカットオフ値より下にある、サンプルはトレーニングセットと大きく異なりません。  図33予測におけるアウトライヤー診断
|
 |
| Software |
 |
 |