確率・統計

(16) 因子分析(Factor Analysis)

「主成分分析(Principal Component Analysis ; PCA)」は、共分散行列が対称行列であることを利用して、「カルーネン・レーベ展開(Karhunen-Loeve Expansion ; KLT)」によってデータを無相関な「主成分」に変換し、その寄与率が高いものを抽出することで、より少ない種類のデータで解析を行うための手法でした。これとよく似た手法として「因子分析(Factor Analysis)」があります。
「因子分析」も、データをより少ない「因子」に変換して分析するという点では「主成分分析」と同様ですが、「主成分分析」がどのようなデータでも(その有効性を無視すれば)適用可能であるのに対し、「因子分析」は統計的モデルが当てはまる場合のみ適用することができます。この章では、「因子分析」について、「主成分分析」との違いも含めて紹介したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 潜在的因子による統計的モデル

ある集団に対して、各要素が複数のデータを持っていたとします。例えば、あるクラスの生徒の、国語・数学・英語の試験の点数などに対し、これをまとめて評価するのはよく行われることです。全教科の和を「総合得点」としたり、数学と語学(国語・英語)の平均から「理系・文系向き」傾向を見るような分析などが考えられます。今、この三教科の得点が一つの因子で説明できると仮定して、次のような統計的モデルを考えます。

zj,i = ajfi + ej,i

j はデータの種別を表し、ここでは三教科 ( j = 1, 2, 3 ) を識別する添字になります。また、i は各生徒に付けられた添字です。zj,i は i 番目の生徒の教科 j の得点で、平均 0、標準偏差 1 に正規化されているものとします。また、fi は各教科に共通な潜在的因子で「共通因子(Common Factor)」と呼ばれ、各生徒に対する fi の値は「因子得点(Factor Score)」といいます。aj は共通因子に対する係数なので、zj,i に対して教科 j がどの程度の重みを持っているかを示す指標になり「因子負荷量(Factor Loading)」、最後の ej,i は、i 番目の生徒の各教科に対する共通因子以外の変動量になるので「独自因子(Unique Factor)」といいます。共通因子と固有因子は、どちらも確率変数であると仮定します。

zj,i の期待値 E[zj,i] は

E[zj,i] = E[ajfi + ej,i] = ajE[fi] + E[ej,i]

と表され、E[zj,i] = 0 より右辺の値もゼロになります。E[fi] と E[ej,i] は平行移動によって右辺がゼロになるように自由に調整することができるので、E[fi] = E[ej,i] = 0 とします。

今、共通因子と独自因子間、また各独自因子間で互いに相関がないものと仮定すると、zj,i の分散 V[zj,i] と共分散 Cov[zj,i,zk,i] ( j ≠ k ) は

V[zj,i]=E[{ ( ajfi + ej,i ) - E[zj,i] }2]
=aj2E[fi2] + 2ajE[fiej,i] + E[ej,i2]
=aj2V[fi] + V[ej,i]
CoV[zj,i,zk,i]=E[{ ( ajfi + ej,i ) - E[zj,i] }{ ( akfi + ek,i ) - E[zk,i] }]
=ajakE[fi2] + ajE[fiek,i] + akE[fiej,i] + E[ej,iek,i]
=ajakV[fi]

と計算することができます。V[zj,i] = 1 より aj2V[fi] + V[ej,i] = 1 であり、V[fi] は aj と V[ej,i] を適当に調整することで 1 とすることができます。また、V[ej,i] = σj2 とすれば、

σj2 = 1 - aj2

で表され、zj,i と zk,i の相関係数 ρj,k

ρj,k = ajak

となり、zj,i の共分散行列 V

V=|1,ρ1,2,ρ1,3|
|ρ2,1,1,ρ2,3|
|ρ3,1,ρ3,2,1|
=|σ12 + a12,a1a2,a1a3|
|a2a1,σ22 + a22,a2a3|
|a3a1,a3a2,σ32 + a32|

で求められます。


因子を複数とした場合、統計的モデルは次のようになります。

zj,i=aj,1f1,i + aj,2f2,i + ... + aj,mfm,i + ej,i
=Σk{1→m}( aj,kfk,i ) + ej,i

上式は、行列とベクトルを使って次のように表すことができます。

|z1,i|=|a1,1,a1,2,...a1,m||f1,i|+|e1,i|
|z2,i||a2,1,a2,2,...a2,m||f2,i||e2,i|
|z3,i||a3,1,a3,2,...a3,m||:||e3,i|
|:||::...:||fm,i||:|
|zp,i||ap,1,ap,2,...ap,m||ep,i|

さらに簡単に

zi = Afi + ei

と表すこともできます。ここでも、共通因子 fk,i と独自因子 ej,i 間、独自因子どうしは互いに無相関であると仮定します。共通因子どうしに対しては二通りの考え方があり、互いに相関がないとした場合の解は「直交解(Orthogonal Solution)」、相関があると仮定した場合は「斜交解(Oblique Solution)」になります。zj,i の分散 V[zj,i] と、zj,i, zl,i の共分散 Cov[zj,i,zl,i] は

V[zj,i]=E[{ [ Σk{1→m}( aj,kfk,i ) + ej,i ] - E[zj,i] }2]
=Σk{1→m}( aj,k2E[fk,i2] ) + 2Σk{1→m-1}( Σl{k+1→m}( aj,kaj,lE[fk,ifl,i] ) ) + 2Σk{1→m}( aj,kE[fk,iej,i] ) + E[ej,i2]
=Σk{1→m}( aj,k2V[fk,i] ) + 2Σk{1→m-1}( Σl{k+1→m}( aj,kaj,lCov[fk,i,fl,i] ) ) + V[ej,i]
Cov[zj,i,zl,i]=E[{ [ Σk{1→m}( aj,kfk,i ) + ej,i ] - E[zj,i] }{ [ Σk{1→m}( al,kfk,i ) + el,i ] - E[zk,i] }]
=Σk1{1→m}( Σk2{1→m}( aj,k1al,k2E[fk1,ifk2,i] ) ) + Σk{1→m}( aj,kE[fk,iel,i] ) + Σk{1→m}( al,kE[fk,iej,i] ) + E[ej,iel,i]
=Σk{1→m}( aj,kal,kV[fk,i] ) + 2Σk1{1→m-1}( Σk2{k1+1→m}( aj,k1al,k2Cov[fk1,i,fk2,i] ) )

であり、特に V[fk,i] = 1 で直交解の場合は

V[zj,i] = Σk{1→m}( aj,k2 ) + V[ej,i]

Cov[zj,i,zl,i] = Σk{1→m}( aj,kal,k )

と非常に単純になります。zj,i が正規化されて分散が 1 ならば、σj2 = V[ej,i] として

Σk{1→m}( aj,k2 ) = 1 - σj2

と表されるので、この値は分散の中で共通因子が占める割合を表すことになり、「共通性(Communality)」と呼ばれる量になります (以下、これを hj2 で表します)。以上の結果から、直交解の場合の共分散行列 V

V=|Σk{1→m}( a1,k2 ) + σ12,Σk{1→m}( a1,ka2,k ),...Σk{1→m}( a1,kap,k )|
|Σk{1→m}( a2,ka1,k ),Σk{1→m}( a2,k2 ) + σ22,...Σk{1→m}( a2,kap,k )|
|::...:|
|Σk{1→m}( ap,ka1,k ),Σk{1→m}( ap,ka2,k ),...Σk{1→m}( ap,k2 ) + σp2|

となり、aj = ( aj,1, aj,2, ... aj,m )T とし、σj2 を対角要素とする行列を D で表せば上式は

V=|a1Ta1,a1Ta2,...a1Tap| + D=|a1T||a1, a2, ... ap| + D = AAT + D
|a2Ta1,a2Ta2,...a2Tap||a2T|
|::...:||:|
|apTa1,apTa2,...apTap||apT|

となります。


主成分分析(Principal Component Analysis ; PCA)」は、「カルーネン-レーベ展開(Karhunen-Loeve Expansion[Transform] ; KLT)」を利用して多変量ベクトルの成分を「主軸変換」によって互いに相関のない「主成分」に変換する手法でした。p 次元の多変量ベクトルを xi ( i = 1, 2, ... N ) としたとき、

V = Σi{1→N}( ( xi - μ )( xi - μ )T ) / N

但し、μ = Σi{1→N}( xi ) / N

は共分散行列になります。最小二乗法より、yi = xi - μ に対して m 個 ( m ≤ p ) の正規直交基底 ej = ( ej,1, ej,2, ... ej,p )T ( j = 1, 2, ... m ) で近似したとき

yi ≅ Σj{1→m}( ci,jej )

但し、ci,j = ( yi, ej )

で表され、さらに、yi 全体に対して最もよく近似できる解を得るためには

Σi{1→N}( ( yi, ej )2 ) = ( ej, Aej )

但し、A = Σi{1→N}( yiyiT )

を最大とするときの ej を求めればよいことになります。これを「カルーネン-レーベ展開」というのでした。A は「半正値対称行列」なので、A の固有値を λj とすれば λj ≥ 0 が成り立ち、λj に対する固有ベクトルを uj = ( u1,j, u2,j, ... up,j ) としたとき Q = ( u1, u2, ... up ) は「直交行列」であり、固有値からなる対角行列を D として固有値分解 A = QDQT が成り立つので

( ej, Aej ) = ( ej, QDQTej ) = ( QTej, DQTej ) = ( e'j, De'j )

但し、e'j = QTej = ( e'j,1, e'j,2, ... e'j,p )T

実際に要素を代入して計算すると

( e'j, De'j )=( e'j,1, e'j,2, ... e'j,p )|λ1,0,...0||e'j,1|
|0,λ2,...0||e'j,2|
|::...:||:|
|0,0,...λp||e'j,p|
=Σk{1→p}( λke'j,k2 )

となります。QT は直交行列なので、ej に対して回転・反転のいずれかの変換を行うのみでノルムは変化させません。( ej, Aej ) は二次形式であり、上式はそれを標準形に変換したことを意味します。共分散行列 V は、上に示した対称行列 A を N で割ったものを意味するので、この手法を使って主軸変換を行うことができます。このとき、固有値分解によって得られる、固有値を対角成分とする対角行列 D は、xi の成分を互いに相関のない状態に変換した場合の共分散行列を意味し、そのときの固有ベクトルが「主成分」となります(*1-1)。

A の固有値分解の式に着目すると、

A= QDQT
=|u1,1,u1,2,...u1,p||λ1,0,...0||u1,1,u2,1,...up,1|
|u2,1,u2,2,...u2,p||0,λ2,...0||u1,2,u2,2,...up,2|
|::...:||::...:||::...:|
|up,1,up,2,...up,p||0,0,...λp||u1,p,u2,p,...up,p|
=|λ1u1,1,λ2u1,2,...λpu1,p||u1,1,u2,1,...up,1|
|λ1u2,1,λ2u2,2,...λpu2,p||u1,2,u2,2,...up,2|
|::...:||::...:|
|λ1up,1,λ2up,2,...λpup,p||u1,p,u2,p,...up,p|
=(λ1u1,λ2u2,...λpup)|u1T|
|u2T|
|:|
|upT|
= Σj{1→p}( λjujujT )

と表すことができます。これは、AujujT による線形結合に変換したことを意味し、その時の係数は λj になります。従って、λj が降順に並べ替えられていて、j が大きくなるほど値が小さくなるとしたとき、m ( < p ) より末尾側の成分を全て切り捨てて

A ≅ Σj{1→m}( λjujujT )

のように A を近似することができて、このとき ( ej, Aej ) は

( ej, Aej ) ≅ Σk{1→m}( λke'j,k2 )

となるので、e'j の m + 1 番目以降の成分は無視されることになります。これはちょうど、p 次元空間内の m 次元超平面上にデータ xi を射影して近似していることを意味します。e'j ( つまり ej ) が m 個あれば、射影したデータは e'j を使って正確に表すことができます。

例として平面上の分布を考えると、二次形式で表現される分布は楕円形となります。楕円形が非常に細い形をしていれば、ある直線上の分布に近似することができます。最も分散が大きくなるような方向の直線に対して射影すれば、この分布は固有値ひとつだけを使って直線上で近似的に表すことができることを意味します。

主成分分析の主軸変換

潜在的因子による統計的モデルでは、複数の因子を共通に表現できる因子を考えます。例えば、二つの因子に対して共通因子を一つ想定したとき、共通因子に比例する形で二つの因子が増減するモデル式を構築します。このとき、i 番目の因子 z1,i, z2,i は、因子負荷量 a1, a2 と共通因子 fi の積に独自因子 e1,i, e2,i を加えた値と等しくなります。これは、回帰式において、fi を独立変数、z1,i, z2,i を従属変数とした場合によく似ていますが、fi も確率変数としていることから、二変数のどちらにも誤差があると仮定した上で、最尤法によって直線への当てはめをした場合に相当します。各都市の人工が fi であれば、z1,i として交通事故の発生件数、z2,i として高齢者が占める割合などが挙げられます。前者は正の相関を持ち、後者は負の相関を持つので、a1 > 0, a2 < 0 であることが予想されるでしょう。E[zj,i] = E[fi] = E[ej,i] = 0 を仮定しているので、回帰式 z = ajf + e は原点を通る直線で、zj,i はその直線付近に分布することになります。zj,i は正規化されているので同じ座標軸上で表現すると、下図のように、原点を通る二つの直線の周囲に zj,i が分布している状態が観察できることになります。

潜在的因子の分布

両者を比較すると、どちらも共分散行列の階数(rank)を近似的に小さくして、ある超平面上に分布を射影する処理を行っていることから同じことをしているように見えます。しかし、主成分分析が、多変量データを互いに相関のないように主軸変換した上である超平面に射影するのに対し、因子分析では、各データを潜在的因子の線形結合で表現できると最初に仮定しています。p 変量のデータに対する主成分分析の場合、各データ xi = ( x1,i, x2,i, ... xp,i ) は主軸変換によって x'i = QTxi に変換され、その第 j 成分 x'j,i

x'j,i = Σk{1→p}( uk,jxk,i )

になります。この逆変換は xi = Qx'i であり、xj,i

xj,i = Σk{1→p}( uj,kx'k,i )

で求められますが、m + 1 ( m < p ) 以降の成分を無視することができれば

xj,i ≅ Σk{1→m}( uj,kx'k,i ) = Σk{1→m}( uj,kΣl{1→p}( ul,jxl,i ) )

と近似することができます。二変量データ ( xi, yi ) で考えれば、

Q=|cosθ,-sinθ|
|sinθ,cosθ|

として

x'i = xicosθ + yisinθ
y'i = -xisinθ + yicosθ

xi = x'icosθ - y'isinθ
yi = x'isinθ + y'icosθ

より、y'i 方向の成分の分散が非常に小さければ

xi ≅ x'icosθ = ( xicosθ + yisinθ )cosθ
yi ≅ x'isinθ = ( xicosθ + yisinθ )sinθ

とすることで、x'i の値だけで xi, yi を表現することができます。

二変量の因子分析では、

xi ≅ a1x'i + e1,i, yi ≅ a2x'i + e2,i

のように共通の因子 x'i で ( xi, yi ) の両方を表すことができると仮定した上でその係数を求めています。主成分分析は、得られたデータの基底を変換して互いに無相関とした上で変換後のデータに意味付けをしているので、統計的モデルを想定しているのではなく結果から意味付けをするのに対し、因子分析は潜在的因子による回帰式が成り立つと仮定した上でその回帰係数を求めていると見ることができるので、最初にモデル式を決めているというところが両者の大きな差異になります。簡単にいえば、主成分は測定データによって構成・説明され、潜在的因子は測定データを構成・意味付けしていることになります。
因子分析は、心理学などで知能の研究に用いられることが多いそうです。知能を測る指標としてはテストの得点を利用するのが一般的であり、テストの得点が潜在的因子となる知能の高さによって説明できるという考え方になります。これを主成分分析に置き換えると、知能の高さがテストの得点によって構成されるという意味になり、本来の考え方と逆転してしまいます。

*1-1)固有値問題(1) 対称行列の固有値」に関する説明を参照


2) 因子の推定

因子分析は統計的モデル式を使って潜在的因子を分析することを目的とした手法であることを説明しました。その具体的な処理内容は次のようになります。

  1. 潜在的因子数の決定
  2. 因子負荷量の推定
  3. 因子軸の回転と因子の解釈
  4. 因子得点の推定

潜在的因子 fk,i の数は、測定データの共分散行列 V の固有値から推測することができます。V は半正値対称行列であり、固有値は必ずゼロ以上であることから、固有値の小さい成分を無視することができる事は前の章で説明したとおりです。目安としては 1.0 をしきい値とすることが多いようです。これは、主成分が一つ以上の成分によって集約されているとの考え方によります。この基準は一般的に「カイザー基準(Kaiser Criterion)」と呼ばれます。その他に、固有値を並べ替えたときの差が非常に大きい箇所で切る方法もあり、これは「スクリー・テスト(Scree Test)」といいます。どの手法を採用するかは実際に固有値を見て判断することになります。

因子数を推定するためのサンプル・プログラムを以下に示します。

/*
  Normalize : データを正規化する

  vector<double>& data : データ列

  戻り値 : True ... 成功 ; False ... 標準偏差がゼロ
*/
bool Normalize( vector<double>& data )
{
  double m = sampleAverage( data );          // 平均値
  double s = sqrt( sampleVariance( data ) ); // 標準偏差

  if ( s == 0 ) {
    cerr << "SD value is zero. ( Data size is zero, or all the data is same. )" << endl;
    return( false );
  }

  for ( unsigned int i = 0 ; i < data.size() ; ++i )
    data[i] = ( data[i] - m ) / s;

  return( true );
}

/*
  createCovarianceMatrix : 共分散行列を求める

  const vector< vector<double> >& f : データ列
  LinearEquationSystem<double>& v : 求めた共分散行列

  戻り値 : True ... 成功 ; False ... データ種の数がゼロ
*/
bool createCovarianceMatrix( const vector< vector<double> >& f, LinearEquationSystem<double>& v )
{
  unsigned int m = f.size();             // データ種の数
  v = LinearEquationSystem<double>( m ); // 共分散行列を表現する行列

  if ( m == 0 ) {
    cerr << "The number of factor is zero." << endl;
    return( false );
  }

  // 共分散行列の要素を計算する
  for ( unsigned int r = 0 ; r < m ; ++r ) {
    v[r][r] = sampleVariance( f[r] );
    for ( unsigned int c = r + 1 ; c < m ; ++c )
      v[r][c] = v[c][r] = sampleCovariance( f[r], f[c] );
  }

  return( true );
}

/*
  calcDiff : データの差分を計算

  const vector<double>& data : 対象データ
  vector<double>& diff : 求めた差分
*/
void calcDiff( const vector<double>& data, vector<double>& diff )
{
  if ( data.size() > 1 ) {
    diff.resize( data.size() - 1 );
    for ( unsigned int i = 0 ; i < diff.size() ; ++i )
      diff[i] = data[i] - data[i + 1];
  } else {
    diff.clear();
  }
}

/*
  printContainer : コンテナクラスの出力

  const vector<double>& data : 対象データ
  vector<double>& diff : 求めた差分
*/
template<class In, class Out> void printContainer( In first, Out last )
{
  if ( first != last ) cout << *( first++ );
  while ( first != last )
    cout << ", " << *( first++ );
  cout << endl;
}

/*
  FA_EstimateFactorCount : 潜在的因子数の推定

  LinearEquationSystem<double> v : 対象の共分散行列
  double qrThreshold : QR法のしきい値
  unsigned int qrMaxCnt : QR法の最大計算回数

  戻り値 : True ... 成功 ; False ... 分散がゼロのデータ列があった, QR法で最大計算回数をオーバーした
*/
bool FA_EstimateFactorCount( LinearEquationSystem<double> v, double qrThreshold, unsigned int qrMaxCnt )
{
  LinearEquationSystem<double> q( 0 ); // 固有値分解後の固有ベクトルからなる直交行列

  cout << "***** Factor Analysis #factor estimation *****" << endl << endl;

  cout << "#class = " << v.size() << endl << endl;

  cout << "***** Covariance Matrix *****" << endl;
  v.printMatrix(); cout << endl;

  // 共分散行列を固有値分解
  Householder_DoubleShiftQR< stat_type >( v, &q, qrThreshold, qrMaxCnt );

  // 固有値分解の結果から固有値を抽出する
  vector<double> eigenValue( v.size() );
  for ( unsigned int i = 0 ; i < v.size() ; ++i )
    eigenValue[i] = v[i][i];
  sort( eigenValue.rbegin(), eigenValue.rend() );

  // 固有値の出力
  cout << "Eigen value\t: ";
  printContainer( eigenValue.begin(), eigenValue.end() );

  // 固有値の差分を計算して出力
  vector<double> d_ev;
  calcDiff( eigenValue, d_ev );
  cout << "Diff\t\t: ";
  printContainer( d_ev.begin(), d_ev.end() );

  // 固有値の二階差分を計算して出力
  vector<double> d2_ev;
  calcDiff( d_ev, d2_ev );
  cout << "2nd Diff\t: ";
  printContainer( d2_ev.begin(), d2_ev.end() );
  cout << endl;

  // スクリー・プロットの出力
  cout << "<<< Scree plot >>>" << endl;
  if ( eigenValue.size() > 0 ) {
    double scale = 5.0;
    for ( double d = eigenValue[0] ; d >= 5.0 ; d /= 5.0 )
      scale /= 5.0;
    for ( unsigned int i = 0 ; i < eigenValue.size() ; ++i )
      cout << i + 1 << "\t:" << string( eigenValue[i] * scale + 1, '*' ) << endl;
  }
  cout << endl;

  return( true );
}

因子分析の実行に先立って、各データ列の平均がゼロ、標準偏差が 1 になるように正規化処理を行うため、関数 Normalize が用意されています。また、データ列から共分散行列を求めるため、関数 createCovarianceMatrix を利用します。
因子数の推定を行う関数が FA_EstimateFactorCount で、自動的に推定した結果を出力するのではなく、求めた固有値やその差分、固有値を降順に並べた上でグラフに表したスクリー・プロット(Scree Plot)などを出力して、因子数の判断は利用者に委ねる形にしています。

共分散行列は対称行列なので、対称行列に特化した固有値計算法 (*2-1) を利用することができます。ここでは「ウィルキンソンの移動法」を利用した QR変換 Householder_DoubleShiftQR を利用しています。

固有値やその差分などは配列に持たせ、その要素を順番に出力するため、関数 printContainer を用意していますが、少しだけ汎用性を持たせて、反復子(Iterator) を持つコンテナクラスなどに対して範囲も含めて指定できるようにしてあります。Standard Template Library(STL) で用意されている反復処理用の関数 for_each と関数オブジェクトを組み合わせた形にする方法もありますが、あまりにも凝りすぎのような気がするのでその方法はあえて避けました。


因子負荷量の推定にはいくつかの方法がありますが、ここでは主因子分析法を紹介します。前にも示した通り、共分散行列 V は、因子負荷量からなる行列 A と独自因子の分散からなる対角行列 D によって次式のように分解することができます。

V = AAT + D

従って、独自因子の分散が分かれば因子負荷量を求めることができますが、通常は D も未知なので、反復法を使って次のような流れで AD を求めることになります。

  1. D の初期値を決める。
  2. 共分散行列 V から D を引いて、その固有値 λj ( j = 1, 2, ... m ) を求める。
  3. 求めた固有値を降順にソートして、因子負荷量の数 m だけ先頭から抽出し、その固有ベクトルを uj とする。
  4. 行列 A = ( √λ1u1, √λ2u2, ... √λmum ) を求める。
  5. 行列 V - AAT を求め、その対角要素からなる行列を D' とする。
  6. DD' を比較して、充分に近ければ処理を終了する。そうでなければ D = D' として 2 から処理を繰り返す。

V - D の m + 1 番目以降の固有値がゼロならば(つまり階数が m ならば)、それは AAT と等しくなります。逆に、求めた AATV から引くと、それは D に等しいことになります。上記処理で求めた A は、m + 1 番目以降の固有値が非常に小さい時に限って D = D' を満たし、処理が完了した段階で AD を得ることができるわけです。

D は対角行列なので対角要素だけ決めれば充分です。その対角要素としては、D の対角要素である σj2 = V[ej,i] の取りうる最大値 1 を使う「ONE 法」や、最小値として知られている、ある種別のデータとそれ以外のデータの間の重相関係数の二乗 ( Squared Multiple Correlation ; SMC ) を使う「SMC法」などがあります ( なお、「(13) 回帰分析法」の「2) 重相関係数(Multiple Correlation Coefficient)」では SMC は「決定係数(Coefficient Of Determination)」または「寄与率(Contribution Ratio)」として紹介しています)。

因子負荷量と独自因子を推定するためのサンプル・プログラムを以下に示します。

/*
  FA_CalcSMC : 各データ列について (分散 - SMC) を算出する(Dの対角成分の初期値)

  *データ列が正規化されていれば 1 - SMC になる。

  const vector< vector<double> >& data : 対象のデータ列
  vector<double>& d : 求めた SMC の配列
*/
void FA_CalcSMC( const vector< vector<double> >& data, vector<double>& d )
{
  unsigned int p = data.size(); // データ種の数

  d.resize( p );
  for ( unsigned int i = 0 ; i < p ; ++i ) {
    const vector<double>& y = data[i]; // i 番目の標本を y とする
    // i 番目の独立変数を除いた標本を作成
    vector< vector<double> > x;
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      if ( j == i ) continue;
      const vector<double>& xj = data[j]; // j 番目の標本
      x.push_back( xj );
    }
    MultipleRegressionCoefficient mrc( x, y );
    d[i] = sampleVariance( y ) - pow( mrc.mcc(), 2 );
  }
}

/*
  FA_EstimateFactorLoading : 因子負荷量の推定(主因子分析法)

  const LinearEquationSystem<double>& v : 共分散行列
  unsigned int m : 因子数
  vector< vector<double> >& a : 因子負荷量
  vector<double>& d : 独自因子(あらかじめ初期値を登録しておく)
  double iterThreshold : 反復法のしきい値
  double qrThreshold : QR法のしきい値
  unsigned int iterMaxCnt : 反復法の最大処理回数
  unsigned int qrMaxCnt : QR法の最大計算回数

  戻り値 : True ... 成功 ; False ... 分散がゼロのデータ列があった, QR/反復法で最大処理回数をオーバーした
*/
bool FA_EstimateFactorLoading( const LinearEquationSystem<double>& v, unsigned int m, vector< vector<double> >& a, vector<double>& d,
                               double iterThreshold, double qrThreshold, unsigned int iterMaxCnt, unsigned int qrMaxCnt )
{
  unsigned int p = v.size(); // データ種の数

  if ( m > p ) {
    cerr << "Specified factor count [m=" << m << "] is bigger than the size of class [p=" << p << "]." << endl;
    return( false );
  }

  if ( p != d.size() ) {
    cerr << "Unigue factor count (" << d.size() << ") is not matched as the class count (" << p << ")." << endl;
    return( false );
  }

  // 因子負荷量のサイズの初期化
  a.resize( p );
  for ( unsigned int i = 0 ; i < p ; ++i )
    a[i].resize( m );

  LinearEquationSystem<double> q( 0 ); // 固有値分解後の固有ベクトルからなる直交行列 Q

  // 反復法の開始
  for ( unsigned int iterCnt = 0 ; iterCnt < iterMaxCnt ; ++iterCnt ) {
    // 共分散行列をコピーして対角成分から独自因子成分を減算 AAt = V - D
    LinearEquationSystem<double> vBuff( v );
    for ( unsigned int i = 0 ; i < p ; ++i )
      vBuff[i][i] -= d[i];

    // 固有値分解 AAt = QD'Q
    Householder_DoubleShiftQR( vBuff, &q, qrThreshold, qrMaxCnt );

    // D'の対角成分から固有値を抽出
    multimap<double,unsigned int> mmapEigenValue; // 固有値と列番号からなるマップ
    for ( unsigned int i = 0 ; i < p ; ++i )
      mmapEigenValue.insert( pair<double,unsigned int>( vBuff[i][i], i ) );

    // 固有値とその固有ベクトルから AAt の対角成分を計算
    vector<double> aat( p, 0 ); // AAt の対角成分
    multimap<double,unsigned int>::const_reverse_iterator cit_ev = mmapEigenValue.rbegin();
    for ( unsigned int i = 0 ; i < m ; ++i ) {
      unsigned int j = cit_ev->second;
      double eigenValue = cit_ev->first;
      for ( unsigned int k = 0 ; k < p ; ++k ) {
        a[k][i] = sqrt( fabs( eigenValue ) ) * q[k][j];
        aat[k] += pow( a[k][i], 2 );
      }
      ++cit_ev;
    }

    // 新たな独自因子成分 D'' を計算し、前回のものと比較
    bool isConverged = true;
    for ( unsigned int i = 0 ; i < vBuff.size() ; ++i ) {
      double new_d = v[i][i] - aat[i];
      if ( pow( d[i] - new_d, 2 ) >= iterThreshold )
        isConverged = false;
      d[i] = new_d; // 新たな D の成分を代入
    }

    if ( isConverged ) return( true );
  }

  cerr << "Run over the max count for iterator method." << endl;

  return( false );
}

V = AAT + D における D の初期値として「SMC法」を使う場合、各対角成分は、対象列のデータの分散から SMC を減算した値になり、データ列が正規化されていれば 1 - SMC で計算できます。この計算を行う関数が FA_CalcSMC です。SMC は前述のとおり、あるデータ種のデータ列とそれ以外のデータ列の間の重相関係数の二乗で求められるので、「2) 重相関係数(Multiple Correlation Coefficient)」で示したサンプル・プログラム内にある MultipleRegressionCoefficient クラスのメンバ関数 mcc が利用できます。

因子負荷量と独自因子を求める関数は FA_EstimateFactorLoading です。ここでは、先に述べた反復法を使い、独自因子が収束するまで処理を繰り返しています。ここでも固有値分解が必要になるため、「ウィルキンソンの移動法」を利用した QR変換 Householder_DoubleShiftQR を利用しています。

ところで、各データ列に対応する独自因子は、D の中の特定の列にあります。ところが、再計算によって得られる D = V - AAT において、AAT は固有値と対応する固有ベクトルを降順に並べて計算しているので、順番が変わってしまうのではないかと感じるかも知れません。
実際には、AAT の値は固有値と固有ベクトルの並べ方には依存しません。A = ( √λ1u1, √λ2u2, ... √λmum ) に対して AAT の値は

AAT = Σj{1→m}( λjujujT )
ujujT =|uj,1||uj,1, uj,2, ... uj,p|=|uj,12,uj,1uj,2,...uj,1uj,p|
|uj,2||uj,2uj,1,uj,22,...uj,2uj,p|
|:||::...:|
|uj,p||uj,puj,1,uj,puj,2,...uj,p2|

で求められるので、r 行 c 列の要素を ar,c とすると

ar,c = Σj{1→m}( λjuj,ruj,c )

となって、固有値・固有ベクトルの位置 j には依存しないことがわかります。また、固有値が小さく無視された固有ベクトルに対しては λjujujT = 0 なので、どの位置にあっても結果には影響しません。計算に必要な値は AAT の対角成分だけなので、サンプル・プログラムの中では対角成分だけを以下の二式で計算しています。

a[k][i] = sqrt( fabs( eigenValue ) ) * q[k][j];
aat[k] += pow( a[k][i], 2 );

a[k][i] は因子負荷量の i 番目の列ベクトルの k 番目の要素を表し、j 番目の固有値 λj と固有ベクトル uj を使って a[k][i] = √λjuj,k と計算しています。q[k][j] は固有ベクトルからなる直交行列の k 行 j 列の要素を表しているので j 列目の固有ベクトルの k 番目の要素という意味になります。k と j が逆になっているように見えますが注意してください。
対角成分については r = c であり、λjuj,ruj,c = λjuj,r2 なので、a[k][i] を二乗することによって AAT の対角成分を求めることができます。


サンプル・プログラムを使って、実際に因子負荷量と独自因子を求めてみたいと思います。以下は、あるクラス 20 人の試験の点数を表した架空のデータです。

表2-1. 20人の五教科テストの点数表
生徒国語数学社会理科英語平均
1588756706968.0
2499450886268.6
3658564978779.6
4517158644257.2
5499060836268.8
6477562817067.0
7416654876362.2
8577058974866.0
9445960644153.6
10407562826364.4
11837157358165.4
12987983718282.6
13793544599863.0
14686049469663.8
15904640576259.0
16708665675167.8
17514556538157.2
18854358817869.0
19815163359565.0
20775853578365.6
平均64.267.357.668.770.765.7

まずは、因子数を推定した結果です。

***** Factor Analysis #factor estimation *****

#class = 5

***** Covariance Matrix *****
( 1, -0.366031, 0.0674933, -0.456982, 0.530843 )
( -0.366031, 1, 0.442199, 0.488902, -0.377569 )
( 0.0674933, 0.442199, 1, 0.197848, -0.0644674 )
( -0.456982, 0.488902, 0.197848, 1, -0.41288 )
( 0.530843, -0.377569, -0.0644674, -0.41288, 1 )

Eigen value     : 2.39129, 1.20284, 0.564695, 0.454689, 0.386478
Diff            : 1.18845, 0.63815, 0.110006, 0.0682103
2nd Diff        : 0.550299, 0.528144, 0.0417955

<<< Scree plot >>>
1       :************
2       :*******
3       :***
4       :***
5       :**

1 以上の値の固有値は二つあり、スクリー・プロットを見てもこの二つが突出していることがわかるので、因子数は 2 として処理を行います。すると次のような結果が得られます。

表2-2. 因子負荷量と独自因子
因子負荷量独自因子共通性
a1a2
国語0.70-0.450.310.69
数学-0.74-0.340.340.66
社会-0.31-0.620.520.48
理科-0.660.00170.570.43
英語0.62-0.210.570.43

この結果をグラフで表したのが下図になります。横軸は一番目の、縦軸は二番目の因子負荷量を示しています。このグラフから、一番目の因子は、語学とそれ以外の教科のどちらが得意であるかを示し、この因子が大きいほど語学の方が得意(逆にそれ以外の教科は苦手)であると解釈できます。また、二番目の因子は全体的な学力を示し、因子負荷量が全体的に負数であることから、この値が小さいほど学力は高いと解釈することができます。

図2-1. 因子負荷量の分布
因子負荷量の分布

因子負荷量の推定には、他にも「セントロイド法(Centroid Method)」や「正準因子分析(Canonical Factor Analysis)」、「最尤法(Maximum-Likelihood Method)」など種々の手法があります。

*2-1)固有値問題 (1) 対称行列の固有値」を参照


3) 因子軸の回転

因子分析で利用される統計的モデル式をもう一度示します。

zi = Afi + ei

ziei はデータ種の数 p 変量のベクトルで、fi は全体的因子の数 m 変量のベクトルです。また、A は因子負荷量を要素とする p 行 m 列の係数行列です。この式は、p 行 p 列の任意の正則行列 T を使って

zi=A(TT-1)fi + ei
=(AT)(T-1fi) + ei
=A*f*i + ei

と表すことができます。但し、A* = AT, f*i = T-1fi となります。この結果は、潜在的因子とその因子負荷量が一意には決まらないことを示していて、因子が解釈しやすい形で変換を行うことができるということを意味します。これを「回転の不定性」といいます。とくに、T が直交行列ならば、その変換は因子軸の回転になるので、解が直交解であれば、変換した結果も直交解のままで、潜在的因子のノルムも変化はしません。

先ほど示した例において、回転行列を使って 30 度回転させた場合の因子は次のようになります。

図・表3-1. 回転後の因子負荷量
因子負荷量
a1a2
国語0.83-0.04
数学-0.47-0.67
社会0.05-0.69
理科-0.57-0.33
英語0.640.13
回転後の因子負荷量の分布

回転後のグラフを見ると、一番目の因子に対しては、語学(国語・英語)と理系(数学・理科)の教科に対する因子負荷量が社会を中心として正負に広がっており、この因子が高いほど語学が得意であると判断することができます。二番目の因子に対しては、語学以外の教科(数学・理科・社会)に対する因子負荷量が負数の値を持っているので、これらの教科が不得意であることを示す因子と解釈することができます。

因子負荷量は、各軸に対して、あるデータ種別のものが絶対値が大きく、それ以外はゼロに近いようにとるのが望ましいとされています。このような形は「単純構造(Simple Structure)」と呼ばれ、先の例では、縦軸に対しては語学の教科、横軸に対しては社会がゼロに近い値となっているので単純構造とみなすことができます。
j 番目の種別の i 番目のデータ zj,i

zj,i = Σk{1→m}( aj,kfk,i ) + ej,i

と表されるのでした。fk,i は i 番目のデータに対する k 番目の潜在的因子、aj,k は j 番目のデータ種別の fk,i ( i = 1, 2, ... ) に対する因子負荷量、ej,i は独自因子をそれぞれ示しています。zj,i は、aj,k のより大きな項に対して強い影響を受けることになるので、単純構造の場合は zj,i の値が特定の潜在的因子のみによって決まり、残りの因子の影響をほとんど受けない状態となることを意味しています。

因子数が二つであれば、グラフを参考に回転を行うことによって最適な解を見つけることができますが、三つ以上になると、各因子軸を二つ選択して回転処理する操作を組み合わせて行う必要があるので、因子数が多くなるほど最適解を得ることが困難になります。そこで、単純構造を機械的に求めるための手法がいくつか提案されています。ここではその中から「バリマックス法(Varimax Method)」を紹介します。

「バリマックス法」では、ある軸に対する因子負荷量の距離の二乗のバラツキができるだけ大きくなるように変換を行います。このとき、因子負荷量の軸からの距離が大きいものと小さいものに分かれ、絶対値の小さい因子負荷量は軸に近い状態になって、単純構造をとることが期待できます。ある軸に対してバラツキを大きくするためには、その分散が大きくなればいいので、

σk2 = Σj{1→p}( aj,k4 ) / p - { Σj{1→p}( aj,k2 ) / p }2

をできるだけ大きくすることを考えます。従って、各軸に対して求めた分散の合計

σ2 = Σk{1→m}( σk2 ) = Σk{1→m}( Σj{1→p}( aj,k4 ) / p - { Σj{1→p}( aj,k2 ) / p }2 )

を最大にするような因子負荷量の変換を行うことによって最適な値を求めることになります。しかし、一度に全ての軸に対して変換を行うことはできないので、実際には m 個の軸の中から二つずつ選んでその中で最大値となるような変換を行う操作を繰り返します。全ての組み合わせに対して処理するためには、mC2 = m( m - 1 ) / 2 回の変換が必要になり、これを一回のパスとして、最大値が収束するまでこのパスを繰り返すことになります。変換処理には、以前「固有値問題 (1) 対称行列の固有値」の「4) ヤコビ法(Jacobi Eigenvalue Algorithm)」で紹介した回転行列による変換処理「ギブンス回転(Givens Rotation)」を利用します。この回転行列は次のようなものでした。

(1)(2)...(t)...(u)...(m)
|1,0,...0,...0,...0,0|(1)
|0,1,...0,...0,...0,0|(2)
|::::::|:
|0,0,...cosθ,...sinθ,...0,0|(t)
|::::::|:
|0,0,...-sinθ,...cosθ,...0,0|(u)
|::::::|:
|0,0,...0,...0,...1,0|
G(i,j,θ) =|0,0,...0,...0,...0,1|(m)

この回転行列を使うことで、aj,t と aj,u ( j = 1, 2, ... p ) を次のように変換することができます。

a'j,t = aj,tcosθ - aj,usinθ

a'j,u = aj,tsinθ + aj,ucosθ

変換後の因子負荷量の分散を σ't2, σ'u2 とすれば、

σ't2 = Σj{1→p}( a'j,t4 ) / p - { Σj{1→p}( a'j,t2 ) / p }2

σ'u2 = Σj{1→p}( a'j,u4 ) / p - { Σj{1→p}( a'j,u2 ) / p }2

この式を θ で微分すると

dσ't2 / dθ = 4Σj{1→p}( a'j,t3( da'j,t / dθ ) ) / p - 2Σj{1→p}( a'j,t2 )[ 2Σj{1→p}( a'j,t( da'j,t / dθ ) ) ] / p2

dσ'u2 / dθ = 4Σj{1→p}( a'j,u3( da'j,u / dθ ) ) / p - 2Σj{1→p}( a'j,u2 )[ 2Σj{1→p}( a'j,u( da'j,u / dθ ) ) ] / p2

da'j,t / dθ と da'j,u / dθ は、

da'j,t / dθ = -aj,tsinθ - aj,ucosθ = -a'j,u

da'j,u / dθ = aj,tcosθ - aj,usinθ = a'j,t

となるので、

dσ't2 / dθ = -4Σj{1→p}( a'j,t3a'j,u ) / p + 4Σj{1→p}( a'j,t2j{1→p}( a'j,ta'j,u ) / p2

dσ'u2 / dθ = 4Σj{1→p}( a'j,u3a'j,t ) / p - 4Σj{1→p}( a'j,u2j{1→p}( a'j,ua'j,t ) / p2

dσ't2 / dθ + dσ'u2 / dθ = 0 のとき、σ't2 + σ'u2 が極値になります。これを Δ とすると、回転によって分散が影響を受けるのは t, u 番目の因子負荷量のバラツキだけでその他は定数なので、Δ = 0 のとき σ2 が極値になることを意味することになって、

Δ=-4Σj{1→p}( a'j,t3a'j,u ) / p + 4Σj{1→p}( a'j,t2j{1→p}( a'j,ta'j,u ) / p2
 + 4Σj{1→p}( a'j,u3a'j,t ) / p - 4Σj{1→p}( a'j,u2j{1→p}( a'j,ua'j,t ) / p2
=-4Σj{1→p}( a'j,t3a'j,u - a'j,u3a'j,t ) / p + 4Σj{1→p}( a'j,t2 - a'j,u2j{1→p}( a'j,ta'j,u ) / p2
=-4Σj{1→p}( a'j,ta'j,u( a'j,t2 - a'j,u2 ) ) / p + 4Σj{1→p}( a'j,t2 - a'j,u2j{1→p}( a'j,ta'j,u ) / p2 = 0

を満たす θ が、σ2 が極値となる回転角になります。

a'j,t2 - a'j,u2=( aj,tcosθ - aj,usinθ )2 - ( aj,tsinθ + aj,ucosθ )2
=aj,t2cos2θ - 2aj,taj,ucosθsinθ + aj,u2sin2θ - aj,t2sin2θ - 2aj,taj,ucosθsinθ - aj,u2cos2θ
=( aj,t2 - aj,u2 )( cos2θ - sin2θ ) - 4aj,taj,ucosθsinθ
=( aj,t2 - aj,u2 )cos2θ - 2aj,taj,usin2θ
a'j,ta'j,u=( aj,tcosθ - aj,usinθ )( aj,tsinθ + aj,ucosθ )
=aj,t2cosθsinθ + aj,taj,ucos2θ - aj,taj,usin2θ - aj,u2cosθsinθ
=aj,taj,u( cos2θ - sin2θ ) + ( aj,t2 - aj,u2 )cosθsinθ
=aj,taj,ucos2θ + ( aj,t2 - aj,u2 )sin2θ / 2
a'j,ta'j,u( a'j,t2 - a'j,u2 )=[ aj,taj,ucos2θ + ( aj,t2 - aj,u2 )sin2θ / 2 ][ ( aj,t2 - aj,u2 )cos2θ - 2aj,taj,usin2θ ]
=aj,taj,u( aj,t2 - aj,u2 )cos22θ - 2aj,t2aj,u2cos2θsin2θ + ( aj,t2 - aj,u2 )2cos2θsin2θ / 2 - aj,taj,u( aj,t2 - aj,u2 )sin2
=aj,taj,u( aj,t2 - aj,u2 )( cos22θ - sin22θ ) + [ ( aj,t2 - aj,u2 )2 / 2 - 2aj,t2aj,u2 ]cos2θsin2θ
=aj,taj,u( aj,t2 - aj,u2 )cos4θ + [ ( aj,t2 - aj,u2 )2 - 4aj,t2aj,u2 ]sin4θ / 4

より

Δ=-4Σj{1→p}( aj,taj,u( aj,t2 - aj,u2 )cos4θ + [ ( aj,t2 - aj,u2 )2 - 4aj,t2aj,u2 ]sin4θ / 4 ) / p
 + 4Σj{1→p}( ( aj,t2 - aj,u2 )cos2θ - 2aj,taj,usin2θ )
 Σj{1→p}( aj,taj,ucos2θ + ( aj,t2 - aj,u2 )sin2θ / 2 ) / p2
=-4cos4θΣj{1→p}( aj,taj,u( aj,t2 - aj,u2 ) ) / p - sin4θΣj{1→p}( ( aj,t2 - aj,u2 )2 - 4aj,t2aj,u2 ) / p
 + 4[ cos2θΣj{1→p}( aj,t2 - aj,u2 ) - sin2θΣj{1→p}( 2aj,taj,u ) ]
 [ cos2θΣj{1→p}( aj,taj,u ) + sin2θΣj{1→p}( aj,t2 - aj,u2 ) / 2 ] / p2

となるので、

A = Σj{1→p}( aj,t2 - aj,u2 )

B = Σj{1→p}( 2aj,taj,u )

C = Σj{1→p}( ( aj,t2 - aj,u2 )2 - 4aj,t2aj,u2 )

D = 4Σj{1→p}( aj,taj,u( aj,t2 - aj,u2 ) )

とすれば式が簡単になって

Δ=-( Dcos4θ + Csin4θ ) / p + 4( Acos2θ - Bsin2θ )( Bcos2θ / 2 + Asin2θ / 2 ) / p2
=-( Dcos4θ + Csin4θ ) / p + 2( ABcos22θ + A2cos2θsin2θ - B2cos2θsin2θ - ABsin22θ ) / p2
=-( Dcos4θ + Csin4θ ) / p + 2[ AB( cos22θ - sin22θ ) + ( A2 - B2 )cos2θsin2θ ] / p2
=-( Dcos4θ + Csin4θ ) / p + 2[ ABcos4θ + ( A2 - B2 )sin4θ / 2 ] / p2

Δ = 0 のとき、cos4θ ≠ 0 ならば、両辺を cos4θ で割れば

-( D + Ctan4θ ) / p + [ 2AB + ( A2 - B2 )tan4θ ] / p2 = 0

より、両辺に p を掛けて辺々整理すると

[ ( A2 - B2 ) / p - C ]tan4θ = D - 2AB / p

よって、

tan4θ = ( D - 2AB / p ) / [ ( A2 - B2 ) / p - C ]

を満たす θ を求めればよいことになります。しかし、極値が最大・最小のいずれとなるかがまだ定まっていないため、Δ をもう一度 θ で微分して

dΔ / dθ=( 4Dsin4θ - 4Ccos4θ ) / p + 2[ -4ABsin4θ + 2( A2 - B2 )cos4θ ] / p2
=4( D - 2AB / p )sin4θ / p - 4[ C - ( A2 - B2 ) / p ]cos4θ / p
=[ 4( D - 2AB / p )sin4θ / p ]{ 1 - [ C - ( A2 - B2 ) / p ] / ( D - 2AB / p )tan4θ }
=4( D - 2AB / p )sin4θ / p{ 1 + [ C - ( A2 - B2 ) / p ]2 / ( D - 2AB / p )2 }

が負になればよいことから

( D - 2AB / p )sin4θ < 0

となるような θ が求める解になります。

求めた因子負荷量をそのまま使って分散の最大値を求める方法は「素バリマックス法(Raw Varimax Method)」と呼ばれます。しかし、この場合は共通性が大きい因子負荷量ほど回転に与える影響が強くなるため、因子負荷量を共通性で割った値を使う手法が考案されています。この方法は「基準バリマックス法Normal Varimax Method」と呼ばれ、一般的にはこの手法が採用されることが多いようです。


バリマックス法による因子の回転を行うためのサンプル・プログラムを以下に示します。

// べき乗を求める関数(doubleのみ)
class Power : public UnaryFunc<double>
{
  double p_;

public:

  Power( double p )
    : p_( p ) {}

  double operator()( const double& x )
  { return( pow( x, p_ ) ); }
};

/*
  transpose : 行列の転置

  const vector< vector<double> >& a : 元の行列
  const vector<double>& h : 共通性の平方根
  vector< vector<double> >& at : 転置した行列

  戻り値 : True ... 転置成功 ; False ... 転置失敗(atは空の行列になる)
*/
bool transpose( const vector< vector<double> >& a, const vector<double>& h, vector< vector<double> >& at )
{
  at.clear();

  unsigned int p = a.size();   // データ種の数
  if ( p < 1 ) return( true ); // a が空なら at も空で正常

  unsigned int m = a[0].size(); // 因子の数
  for ( unsigned int i = 1 ; i < p ; ++i ) {
    if ( a[i].size() != m ) {
      cerr << i << "th column size of matrix is not matched." << endl;
      return( false );
    }
  }

  // 転置処理開始(共通性での除算を実施)
  at.resize( m );
  for ( unsigned int i = 0 ; i < m ; ++i ) {
    at[i].resize( p );
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      if ( h[j] != 0 ) at[i][j] = a[j][i] / h[j];
    }
  }

  return( true );
}

/*
  calcMaxVarAngle : 分散を最大にする回転角を求める

  vector<double>& a0, vector<double>& a1 : データ列

  戻り値 : 回転角(i0,i1がデータ列のサイズ以上なら NAN を返す)
*/
double calcMaxVarAngle( vector< vector<double> >& at, unsigned int i0, unsigned int i1 )
{
  if ( at.size() == 0 ) return( NAN );
  if ( i0 >= at.size() || i1 >= at.size() ) return( NAN );

  unsigned int p = at[0].size(); // データ種のサイズ

  // 回転角を求めるための項(A,B,C,D)を求めるための計算
  vector<double> vec_a( p ), vec_b( p ), vec_c( p ), vec_d( p );
  for ( unsigned int i = 0 ; i < p ; ++i ) {
    vec_a[i] = ( at[i0][i] + at[i1][i] ) * ( at[i0][i] - at[i1][i] );
    vec_b[i] = at[i0][i] * at[i1][i] * 2;
    vec_c[i] = ( vec_a[i] + vec_b[i] ) * ( vec_a[i] - vec_b[i] );
    vec_d[i] = vec_a[i] * vec_b[i] * 2;
  }

  // 回転角を求めるための項(A,B,C,D)
  double asum = sum( vec_a );
  double bsum = sum( vec_b );
  double csum = sum( vec_c );
  double dsum = sum( vec_d );

  double numerator = dsum - 2 * asum * bsum / (double)p; // D - 2AB / p
  double theta = atan( numerator / ( ( asum + bsum ) * ( asum - bsum ) / (double)p - csum ) ); // 回転角(4θ)

  // 符号が逆の場合、π 加える
  if ( numerator * sin( theta ) > 0 )
    theta += M_PI;

  return( theta / 4.0 );
}

/*
  FA_Varimax : 因子負荷量の回転(基準Varimax法)

  vector< vector<double> >& a : 因子負荷量
  double threshold : 反復法のしきい値
  unsigned int maxCnt : 反復法の最大処理回数

  戻り値 : True ... 成功 ; False ... 反復法で最大処理回数をオーバーした
*/
bool FA_Varimax( vector< vector<double> >& a, double threshold, unsigned int maxCnt )
{
  unsigned int p = a.size(); // データ種の数
  if ( p < 1 ) return( true );
  unsigned int m = a[0].size(); // 潜在的因子の数

  Power pow4( 4 ); // 四乗和計算用関数オブジェクト
  Power pow2( 2 ); // 二乗和計算用関数オブジェクト

  vector<double> h( p ); // 共通性の平方根
  for ( unsigned int i = 0 ; i < p ; ++i )
    h[i] = sqrt( sum( a[i], pow2 ) );

  vector< vector<double> > at; // a の二乗の転置行列
  if ( ! transpose( a, h, at ) ) {
    cerr << "Failed to transpose factor loading matrix." << endl;
    return( false );
  }

  vector<double> vk( m ); // 各列の分散
  for ( unsigned int i = 0 ; i < m ; ++i )
    vk[i] = sum( at[i], pow4 ) / (double)p + pow( sum( at[i], pow2 ) / (double)p, 2 );
  double v = sum( vk ); // 分散の合計

  // 反復処理開始
  for ( unsigned int cnt = 0 ; cnt < maxCnt ; ++cnt ) {
    for ( unsigned int i = 0 ; i < m ; ++i ) {
      for ( unsigned int j = i + 1 ; j < m ; ++j ) {
        double theta = calcMaxVarAngle( at, i, j ); // 分散を最大にする回転角
        for ( unsigned int k = 0 ; k < p ; ++k ) {
          double ai = a[k][i] * cos( theta ) - a[k][j] * sin( theta ); // k 行 i 列の更新
          double aj = a[k][i] * sin( theta ) + a[k][j] * cos( theta ); // k 行 j 列の更新
          a[k][i] = ai;
          a[k][j] = aj;
          if ( h[k] != 0 ) {
            at[i][k] = ai / h[k];
            at[j][k] = aj / h[k];
          }
        }
      }
    }
    for ( unsigned int i = 0 ; i < m ; ++i )
      vk[i] = sum( at[i], pow4 ) / (double)p + pow( sum( at[i], pow2 ) / (double)p, 2 );
    double new_v = sum( vk ); // 新たに求めた分散

    if ( fabs( v - new_v ) < threshold )
      return( true );

    v = new_v;
  }

  cerr << "Run over the max count for iterator method." << endl;
  return( false );
}

バリマックス法による因子の回転は FA_Varimax で行っています。ここでは共通性で因子負荷量を割った上で分散を求める「基準バリマックス法」を用いています。FA_Varimax の中で、因子負荷量の行列 a の要素を二乗した値を持つ転置行列 at を作成する処理を行っています。分散を計算するときは列ベクトル単位(因子ごと)であり、a は行ベクトル単位(データ種ごと)に保持しているため、計算を簡単に行えるようにするためです。転置行列を作成する処理は関数 transpose で行い、さらに共通性で除算する処理もまとめて行っています。

分散を最大にする回転角を求める処理は calcMaxVarAngle で行っています。処理の内容は先に説明した内容そのままですが、sinθ の符号を逆転するときに θ に π を加算している部分について若干補足しておくと、ここでは tanθ の符号を変えずに sinθ の符号を逆転させる必要があり、π を加算すると、

sin( θ + π ) = -sinθ

cos( θ + π ) = -cosθ

より tan( θ + π ) = tanθ なので、tanθ の符号を変えずに sinθ の符号のみを逆転させることができます。しかし -1 を掛けた場合は cos( -θ ) = cosθ なので tanθ の符号も逆転してしまい、この方法は利用できません。


バリマックス法を使い、先に示した五教科の試験結果に対する因子負荷量を処理すると、回転角は 28.9 度になり以下のような結果が得られます。これは、先ほどグラフを参考にして処理したときの回転角と非常に近い値となっています。

表3-2. バリマックス法により回転した因子負荷量
a1a2
国語0.83-0.055
数学-0.48-0.66
社会0.032-0.70
理科-0.58-0.32
英語0.650.11

バリマックス法は「ヘンリー・カイザー(Henry Felix Kaiser)」によって 1958 年に提案された手法で、直交行列による変換なので、因子負荷量が直交解の場合は変換後もやはり直交解のままとなります。直交変換を行う手法としては他に「コーティマックス法(Quartimax Method)」や「エカマックス法(Equamax Method)」などがあります。
直交軸への変換ではなく、斜交軸へ変換する手法もあり、「オブリマックス法(Oblimax Method)」や「オブリミン法(Oblimin Method)」などが有名です。これらは斜交解へ変換することになり、各因子に対して相関があるような解が得られます。


4) 因子得点の推定

因子負荷量 aj,k が推定できたので、データ zj,i

zj,i = Σk{1→m}( aj,kfk,i ) + ej,i

で表されることを利用して因子得点 fk,i が決まります。独自因子 ej,i は各共通因子には依存せず、平均ゼロの確率変数であることから、いわば誤差因子と見ることができるので、ej,i がゼロの時の因子得点の推定値を f~k,i としたとき、zj,i を j = 1 から p まで並べた連立方程式は、行列を使って

zi = Af~i

と表されます。ここで、zi は p 次元の列ベクトル、A は p 行 m 列の行列 ( 但し p ≥ m )、f~i は m 次元の列ベクトルとなります。p ≥ m のとき、A の階数(rank)が m ならば、LA = E ( 但し、E は m 行 m 列の単位行列 ) を満たす m 行 p 列の行列 L (これを左逆行列といいます)が存在するので(*4-1)、両辺に L を掛けることで

Lzi = LAf~i = f~i

L の k 行 j 列目の要素を bk,j とすれば、k 行目の因子得点 f~k,i

f~k,i = Σj{1→p}( bk,jzj,i )

になります。k 番目のデータ種に対し、bk,j は p 個あり、それに対して連立方程式の数は、データ数を N としたとき N 個です。通常は p < N であることから、bk,j を未知数としたとき、未知数に対して連立方程式の数の方が多く、一般的には解は不定になります。従って、最小二乗法を利用して、bk,j の最尤値を求めることになります。

fk,i と f~k,i との差の二乗和の 1 / 2 を Jk とすると、

Jk=Σi{1→N}( ( fk,i - f~k,i )2 ) / 2
=Σi{1→N}( ( fk,i - Σj{1→p}( bk,jzj,i ) )2 ) / 2

∂Jk / ∂bk,s

∂Jk / ∂bk,s=Σi{1→N}( -zs,i[ fk,i - Σj{1→p}( bk,jzj,i ) ] )
=i{1→N}( zs,ifk,i ) + Σi{1→N}( Σj{1→p}( bk,jzs,izj,i ) )
=i{1→N}( zs,ifk,i ) + Σj{1→p}( bk,jΣi{1→N}( zs,izj,i ) )

よって、∂Jk / ∂bk,s = 0 のとき

Σj{1→p}( bk,jΣi{1→N}( zs,izj,i ) ) = Σi{1→N}( zs,ifk,i )

が成り立ちます。左辺にある

Σi{1→N}( zs,izj,i )

zi の s 番目と j 番目の要素に対する標本共分散の N 倍であり、zi の各要素が平均ゼロ、標準偏差 1 に正規化されていれば、これは相関係数の N 倍を表すことになるので、s 番目と j 番目の要素に対する相関係数を rs,j で表せば

(左辺) = Σj{1→p}( bk,j・Nrs,j )

となります。右辺は

(右辺)=Σi{1→N}( zs,ifk,i )
=Σi{1→N}( fk,i[ Σt{1→m}( as,tft,i ) + es,i ] )
=Σi{1→N}( fk,iΣt{1→m}( as,tft,i ) ) + Σi{1→N}( fk,ies,i )
=Σt{1→m}( as,tΣi{1→N}( fk,ift,i ) ) + Σi{1→N}( fk,ies,i )

と変形できて、第一項の

Σi{1→N}( fk,ift,i )

は、因子得点が fk,i が互いに無相関(直交解)で、分散が 1 に正規化されていれば Nδkt ( δkt : クロネッカのデルタ ; = 0 ( k ≠ t ), = 1 ( k = t ) )で表すことができます。また、第二項は、共通因子と独自因子が互いに無相関であるという仮定からほぼゼロであるとみなすことができて、

(右辺) = Σt{1→m}( as,t・Nδkt ) = Nas,k

となるので、正規方程式は

Σj{1→p}( rs,jbk,j ) = as,k ( s = 1, 2, ... p ; k = 1, 2, ... m )

と求められ、行列で表せば

|r1,1,r1,2,......r1,p||b1,1b2,1...bm,1|=|a1,1,a1,2,...a1,m|
|r2,1,r2,2,......r2,p||b1,2b2,2...bm,2||a2,1,a2,2,...a2,m|
|::......:||::...:||::...:|
|::......:||::...:||::...:|
|rp,1,rp,2,......rp,p||b1,pb2,p...bm,p||ap,1,ap,2,...ap,m|

となるので、相関係数 rs,j を並べた相関行列を R、因子負荷量からなる行列を A、未知数 bs,k からなる行列を B とすれば、正規方程式は

BT = R-1A

であり、R-1 の成分を rs,j で表せば

bk,s = Σj{1→p}( rs,jaj,k ) ( s = 1, 2, ... p ; k = 1, 2, ... m )

より因子得点の推定値は

f~k,i=Σs{1→p}( bk,szs,i )
=Σs{1→p}( Σj{1→p}( rs,jaj,k )zs,i )

になります。


実際のデータを使って因子得点の推定値を求めてみます。下表は、前に示した 20 人分のテストの得点表を正規化した値を表しています。

表4-1. 20人の五教科テストの点数表(正規化後)
生徒国語数学社会理科英語平均
1-0.351.17-0.190.072-0.100.12
2-0.871.59-0.881.07-0.510.080
30.0491.050.741.580.960.88
4-0.760.220.046-0.26-1.69-0.49
5-0.871.350.280.80-0.510.21
6-0.990.460.510.69-0.0410.13
7-1.33-0.077-0.421.02-0.45-0.25
8-0.410.160.0461.58-1.340.0065
9-1.16-0.490.28-0.26-1.75-0.68
10-1.390.460.510.74-0.45-0.027
111.090.22-0.070-1.880.61-0.0066
121.950.702.950.130.671.28
130.85-1.93-1.58-0.541.61-0.32
140.22-0.44-1.00-1.261.49-0.20
151.49-1.27-2.04-0.65-0.51-0.60
160.341.110.86-0.095-1.160.21
17-0.76-1.33-0.19-0.870.61-0.51
181.20-1.450.0460.690.430.18
190.97-0.970.63-1.881.430.036
200.74-0.55-0.53-0.650.73-0.055

このデータの共分散行列 R は次のようになります。

R =|1,-0.37,0.067,-0.46,0.53|
|-0.37,1,0.44,0.49,-0.38|
|0.067,0.44,1,0.20,-0.064|
|-0.46,0.49,0.20,1,-0.41|
|0.53,-0.38,-0.064,-0.41,1|

また、R の逆行列は R-1 は次のようになります。

R-1 =|1.67,0.35,-0.39,0.42,-0.61|
|0.35,1.75,-0.69,-0.46,0.24|
|-0.39,-0.69,1.35,-0.11,-0.014|
|0.42,-0.46,-0.11,1.53,0.23|
|-0.61,0.24,-0.014,0.23,1.51|

バリマックス回転後の因子負荷量は次のような結果でした。

表4-2. 因子負荷量(回転後)
aj,1aj,2
国語 a1,k0.83-0.055
数学 a2,k-0.48-0.66
社会 a3,k0.032-0.70
理科 a4,k-0.58-0.32
英語 a5,k0.650.11

これらの結果から、未知数 bk,s を要素とする m x p 行列 B の 1 行 1 列目の要素 b1,1 は次のような式で求めることができます。

b1,1=Σj{1→p}( r1,jaj,1 )
=1.67 x 0.83 + 0.35 x (-0.48) + (-0.39) x 0.032 + 0.42 x (-0.58) + (-0.61) x 0.65
=0.57

B の全要素 bk,s の計算結果を以下に示します。

表4-3. bk,s の計算結果
bk,1bk,2bk,3bk,4bk,5
b1,s0.57-0.150.11-0.180.23
b2,s-0.25-0.52-0.44-0.11-0.024

B から因子得点の推定量 f~k,i を求めることができます。例えば、1 番目の生徒の第一因子得点の推定値 f~1,1 は次のように計算できます。

f~1,1=Σs{1→p}( b1,szs,1 )
=0.57 x -0.35 + (-0.15) x 1.17 + 0.11 x (-0.19) + (-0.18) x 0.072 + 0.23 x (-0.10)
=-0.43

因子得点の推定量 f~k,i を求めた結果を以下に示します。

表4-4. 因子得点の推定量
f~1,if~2,i
f~k,1-0.43-0.44
f~k,2-1.13-0.32
f~k,3-0.11-1.07
f~k,4-0.800.12
f~k,5-0.92-0.67
f~k,6-0.70-0.28
f~k,7-1.070.46
f~k,8-0.83-0.13
f~k,9-0.910.49
f~k,10-1.03-0.18
f~k,111.04-0.17
f~k,121.45-2.16
f~k,131.061.49
f~k,140.640.70
f~k,150.811.25
f~k,16-0.13-1.00
f~k,170.0401.03
f~k,180.880.34
f~k,191.420.15
f~k,200.720.38

一番目の因子は語学が得意かどうかを示し、二番目の因子は語学以外の分野が苦手であるかどうかを示していると解釈できるのでした。1 から 10 番目までの前半の生徒の、一番目の因子得点 f~1,i は、後半の生徒に比べて低い値(負値)を示しており、これは、国語と英語の点数が相対的に低いという事実と一致しています。また、12 番目の生徒は f~1,12 の値が大きく、逆に f~2,12 の値は小さいので、全教科について得意であると解釈ができて、実際のテストの点数を見ても、合計点数は最も高く、バラツキもそれほどないことがわかります。


因子得点の推定量を求めるためのサンプル・プログラムを以下に示します。

/*
  FA_EstimateFScoreCoef : 因子得点の線形式の係数を求める

  LinearEquationSystem<double>& v : 共分散行列
  const vector< vector<double> >& a : 因子負荷量
  vector< vector<double> >& b : 求めた係数

  戻り値 : True ... 成功 ; False ... 反復法で最大処理回数をオーバーした
*/
bool FA_EstimateFScoreCoef( LinearEquationSystem<double>& v, const vector< vector<double> >& a, vector< vector<double> >& b )
{
  unsigned int p = v.size(); // データ種のサイズ(vのサイズ,aの行数)
  if ( p < 1 ) {
    cerr << "The size of covariance matrix seems to be zero." << endl;
    return( false );
  }
  if ( a.size() != p ) {
    cerr << "Row size of Factor Loading(" << a.size() << ") is not same as the size of covariance matrix(" << p << ")." << endl;
    return( false );
  }

  unsigned int m = a[0].size(); // 因子負荷量の数
  for ( unsigned int i = 1 ; i < p ; ++i ) {
    if ( a[i].size() != m ) {
      cerr << "Column size of " << i << "th row(" << a[i].size() << ") is not same as other row(" << m << ")." << endl;
      return( false );
    }
  }

  LinearEquationSystem<double> inv( p ); // vの逆行列
  if ( ! Inverse( v, inv ) ) {
    cerr << "Failed to determine an inverse of covariance matrix." << endl;
    return( false );
  }

  b.resize( m );
  for ( unsigned int i = 0 ; i < m ; ++i ) {
    b[i].resize( p );
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      vector<double> buff( p ); // 和を計算するためのバッファ
      for ( unsigned int k = 0 ; k < p ; ++k ) {
        buff[k] = inv[j][k] * a[k][i];
      }
      b[i][j] = sum( buff );
    }
  }

  return( true );
}

/*
  FA_EstimateFactorScore : 因子得点の推定値を求める

  const vector<double>& b : 因子得点の線形式の係数
  const vector<double>& d : データ列

  戻り値 : 因子得点の推定値
*/
double FA_EstimateFactorScore( const vector<double>& b, const vector<double>& d )
{
  unsigned int p = b.size(); // データ種の数
  if ( p != d.size() ) {
    cerr << "The size of data(" << d.size() << ") is not same as the size of coefficient(" << p << ")." << endl;
    return( NAN );
  }

  vector<double> buff( p ); // 和を計算するためのバッファ
  for ( unsigned int i = 0 ; i < p ; ++i )
    buff[i] = b[i] * d[i];

  return( sum( buff ) );
}

/*
  FA_PrintFactorScore : 各データの因子得点の推定値を計算して出力する

  const vector< vector<double> >& b : 因子得点の線形式の係数(m x p 行列)
  const vector<double>& f : データ列(p x n 行列)

  m : 因子の数 ; p : データ種の数 ; n : データ数
*/
void FA_PrintFactorScore( const vector< vector<double> >& b, const vector< vector<double> >& f )
{
  unsigned int p = f.size();    // データ種の数
  if ( p == 0 ) return;

  unsigned int m = b.size();    // 因子負荷量の数
  unsigned int n = f[0].size(); // データ数

  vector<double> z( p ); // f の列を抽出したベクトル
  for ( unsigned int i = 0 ; i < n ; ++i ) {
    for ( unsigned int r = 0 ; r < p ; ++r )
      z[r] = f[r][i];
    for ( unsigned int j = 0 ; j < m ; ++j ) {
      if ( j > 0 ) cout << ",";
      cout << FA_EstimateFactorScore( b[j], z );
    }
    cout << endl;
  }
}

未知数である係数 B を求めるための関数が FA_EstimateFScoreCoef で、共分散行列 v と因子負荷量 a を渡すことで、係数を b に返します。ここで求めた係数 b を元に、FA_EstimateFactorScore を使って因子得点の推定量を求めることができます。この関数で一度に求めることができるのは、集団の中のある要素が持つデータ列に対する一つの因子得点のみで、これを各要素、因子ごとに繰り返し行う必要があります。この処理を繰り返して、結果を出力するための関数が FA_PrintFactorScore です。

*4-1)確率・統計 (12) 二標本の解析 -2-」の「補足2) べき等行列の階数」に関する説明を参照


因子分析については、因子負荷量の推定法や因子の回転・解釈の方法にいく通りものやり方があり、今回は、その中で直交解を求めるための代表的な手法だけを紹介しました。また機会があれば、斜交解を求めるための手法なども紹介していきたいと思います。


<参考文献>
1. 「多変量統計解析法」 田中 豊・脇本和晶 共著 (現代数学社)
2. 愛知学院大学
「千野研究室」様のホームページ - SASによるデータ解析/基礎と応用
3. Wikipedia

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system