先日の真夜中、データを図示するという地味な作業をコリコリとやっていて、あれ？ これってなんだっけ？ と混乱してしまったことがあった。
個体変量行列でも集計表でもなんでもいいけど、とにかく行列\(X\)があるとするじゃないですか。これを列中心化して、主成分分析やって、各変数(\(X\)の列)の第一主成分負荷と第二主成分負荷を座標にして、変数を二次元空間にマッピングしたとしますよね。同時に、各ケース(\(X\)の行)の第一主成分得点と第二主成分得点を座標にして、ケースを二次元空間にマッピングしたとしますよね。
そこに誰かがやってきて、ちょっとちょっと、\(X\)と同じ列を持ってる追加のケースがあるんだけどさ、その図に無理やり乗っけてくれない？ と云ってきたとしよう。いいですよってんで、主成分負荷を変えずにその追加ケースの主成分得点を出してケースの図に乗せてあげる。こういうことはよくある。これをsupplementary caseなどという。
また別の誰かがやってきて、ちょっとちょっと、\(X\)と同じ行を持ってる追加の変数があるんだけどさ、その図に無理やり乗っけてくれない？ と云ってきたとしよう。いいですよってんで、この変数の主成分負荷を出して変数の図に乗せてあげる。こういうこともよくある。これをsupplementary variableなどという。
ここで、あれれ、と混乱してしまったのである。俺なにやってんだ？
いや、やり方がわからないという話ではない。Rであれば、追加ケースの得点はprcomp()が返すオブジェクトに対してpredict()すれば出てくる。追加変数の負荷はprcomp()では出せないと思うけど、たとえばfactominer::PCA()で出せる。
理屈が分からないという話でもない。列中心化済行列\(X\)のSVDが\(X = U D V^\top\)だとして、ここからなんらかの形で、\(FG^\top=X\)となるように、ケースの座標行列\(F\)と変数の座標行列\(G\)を作ったとする (たとえば\(F=U, G=VD\)。いわゆるGabrielの共分散biplotがそうだ)。追加ケースの座標とは、追加ケースから\(X\)の列平均を引いた行列\(X_{suprow}\) について \(X_{suprow} = F_{suprow} G^\top\)となる\(F_{suprow}\)だろうし、追加変数の座標とは、追加変数を列中心化した行列\(X_{supcol}\)について\(X_{supcol} = F G_{supcol}^\top\)となる\(G_{supcol}\)であろう。
そういう話をしたいわけじゃなくて… 追加変数の負荷ってなんなの？ いったいどういう意味を持っているの？ と思ったのである。追加ケースの得点はわかる。それはもとの主成分分析とまったく同じ重みを使って求めた合成得点だ。では追加変数の負荷ってなに？ 感覚としては、もとの変数の数をp1, 追加変数の数をp2として、p1+p2次元の散布図をぐるぐる回してその影を障子に映し(嗚呼、なんと文学的な表現だろうか)、もとのp1個の変数について影が広がって見えるようないい感じの角度で止めたとき、p2本の軸がなす障子紙に対してなす角度、みたいなもんであろう。オーケー、でもそれは結局なにを表しているの？ 俺はこんな真夜中にいったい何をやっているの？
Graffelman, J., Aluja-Banet, T. (2003) Optimal Representation of Supplementary Variables in Biplots from Principal Component Analysis and Correspondence Analysis. Biometrical Journal, 45(4), 491-509.
… その晩は頭に疑問符を抱えながら、とにかくチャートを描いて送って、別の仕事で徹夜して、翌日昼にふらふらになって寝た。それから数日ぼーっと過ごし、ふと思い出して、それらしい題名の解説論文を探して読んでみた。
著者らはバルセロナの大学の先生。google様いわく被引用件数31件。
1. イントロダクション
主成分分析(PCA)や対応分析(CA)では、元の分析に用いなかったケースを補助点(supplementary points)として表示することがある。求め方はGreenacre(1984, “Theory and Applications of CA”)のp.70, Greenacre(1993, “CA in Practice”)のChap.12をみよ。元の分析に用いなかった変数を補助変数(supplementary variables)として表示することもある。環境研究ではこれをindect gradient analysisという[←へー]。
まあとにかく、どちらも基本的な公式はGabriel(1995 Chap.)で与えられていた。本論文では、それらの公式を幾何的に正当化する。つまり、その最小二乗的な意味での最適性を示す。
2. 補助変数と補助点
もとのデータ行列を\(\mathbf{X}\) (サイズ\(n \times p\), ランク\(r\))とする。SVDをつかって\(\mathbf{X}=\mathbf{FG’}\)となんらか分解して、\(\mathbf{F}\)と\(\mathbf{G}\)のそれぞれ最初の2列を座標にしてバイプロットを描きました、と。
ある補助変数\( \mathbf{z} \)(サイズ\(n \times 1\))をこのバイプロット上のベクトル\(\mathbf{v}\)で表すことを考える。ケース\(i\)の真の値を\(z_i\), バイプロットから復元した値を\(\hat{z}_i\)として、\(\mathbf{e} = \mathbf{z} – \mathbf{\hat{z}}\)の二乗和を最小化したい。
つまり、\(\mathbf{\hat{z}} = \mathbf{Fv}\)と近似して、$$ \mathbf{e’ e} = (\mathbf{z} – \mathbf{Fv})’ (\mathbf{z} – \mathbf{Fv})$$を最小化したいわけですよね。その解は$$\mathbf{v} = (\mathbf{F’ F})^{-1} \mathbf{F’z}$$ つまり \(\mathbf{z}\)を\(\mathbf{F}\)の列に回帰したときの回帰係数である。[あっ、そうか、なるほど]
こんどは、あるケース\(\mathbf{x}\) (サイズ\(p \times 1\))をこのバイプロット上の補助点\(\mathbf{p}\)として表すことを考える。$$ \mathbf{e’ e} = (\mathbf{x} – \mathbf{G p})’ (\mathbf{x} – \mathbf{G p}) $$を最小化する\(p\)は $$ \mathbf{p} = (\mathbf{G’G})^{-1} \mathbf{G’ x}$$ つまり\(\mathbf{x}\)を\( \mathbf{G}\)の列に回帰したときの回帰係数である。
一般化しよう。補助情報\(Z\) (サイズ\(m \times q\))の「条件付きバイプロット」を構築するという問題を考える。条件付きバイプロットとは、行だか列だかの座標が固定されているバイプロットである。
\(\mathbf{X} = \mathbf{F G’}\)としよう。補助変数問題とは\(m=n, \mathbf{A} = \mathbf{F}\)として\(\mathbf{Z} = \mathbf{A B’}\)となる\(\mathbf{B}\)を探せという問題であり、解は$$ \mathbf{B’} = (\mathbf{F’F})^{-1} \mathbf{F’Z} $$ である。補助点問題とは\(q=p, \mathbf{B} = \mathbf{G}\)として\(\mathbf{Z} = \mathbf{A B’}\)となる\(\mathbf{A}\)を探せという問題であり、解は$$ \mathbf{A’} = (\mathbf{G’G})^{-1} \mathbf{G’Z’} $$ である。
3. 主成分分析
\(\mathbf{X}\)は中心化済みだとして、$$ \frac{1}{\sqrt{n}} \mathbf{X} = \mathbf{U D V’}$$とSVDしますよね。\(\sqrt{n}\)で割ったのは特異値が\(\mathbf{X}\)の共分散行列の固有値の平方根と等しくなるようにしたかっただけです。
バイプロットの行座標\( \mathbf{F} \)と列座標\( \mathbf{G} \)の決め方はいろいろありえます。二つのやり方をご紹介しましょう。$$ \mathbf{F}_s = \sqrt{n} \mathbf{U}, \ \ \mathbf{G}_p = \mathbf{V D} $$ $$ \mathbf{F}_p = \sqrt{n} \mathbf{U D}, \ \ \mathbf{G}_s = \mathbf{V} $$ \(p\)はprincipal coordinatesのp、\(s\)はstandard coordinatesのsです。これはGreenacreさんの呼び方です。主成分分析の場合でいうと、\(\mathbf{F}_p\)にはいっているのは「主」成分で、\(\mathbf{F}_s\)にはいっているのは「標準化された」主成分なわけです。覚えやすいですね。もっとも\(\mathbf{G}_s\)にはいっている列座標は全然標準化されてないですけどね。混乱しますね。
[ほんとにね、主座標と標準座標ってどっちがどっちだかいっつもわかんなくなる。特異値\(\mathbf{D}\)を掛けるほうが主座標で掛けないほうが標準座標ね、よし覚えた]
\(\mathbf{X}\)が標準化済みであれば、\(\mathbf{G}_p\)は主成分と変数との相関行列である。\(\mathbf{F}_s\)と\(\mathbf{G}_p\)のバイプロットは、ベクトルのなす角のコサインが相関を近似し、ベクトルの長さは主成分への回帰によって説明された分散の平方根を表す。いっぽう\(\mathbf{F}_p\)と\(\mathbf{G}_s\)のバイプロットでは、行の点のユークリッド距離がケース間のユークリッド距離を近似する。
さて、補助変数について考えよう。
\(\mathbf{F}_s\)と\(\mathbf{G}_p\)のバイプロットなら $$ \mathbf{v} = \frac{1}{n} \mathbf{F’}_s \mathbf{z} $$ となるわけね。\(\mathbf{z}\)が標準化済みなら、これは主成分と補助変数の相関係数ベクトルである。バイプロットで説明された分散は[…導出を省略…] 主成分との相関の二乗和 \( \sum_k r^2_k(\mathbf{z}, \mathbf{F}_{sk}\)である。二次元のバイプロットなら、補助変数のベクトルの長さは、その補助変数を2つの主成分に回帰したときの決定係数の平方根になる。つまり補助変数のベクトルは、分析に使った変数のベクトルと同じように、単位円の内側に落ちるわけである。
\(\mathbf{F}_p\)と\(\mathbf{G}_s\)のバイプロットならどうなるか。このとき$$ \mathbf{v} = \mathbf{D}^{-1} \left( \frac{1}{n} \mathbf{F’}_s \mathbf{z} \right) $$ である。\(\mathbf{v}\)は主成分と補助変数の相関係数ベクトルではなく、それを特異値で割ったものである。補助変数のベクトルの長さにはもはやあまり意味がない。
\(k\)水準のカテゴリカルな補助変数を\(k\)個の指標変数にコーディングしたとき、\(\mathbf{F}_s\)と\(\mathbf{G}_p\)のバイプロットなら、それらのベクトルはそのカテゴリについて2つの標準化された主成分を平均したものになる。\(\mathbf{F}_p\)と\(\mathbf{G}_s\)のバイプロットではそうならない。
ついでに補助点についても考えておくと、\(\mathbf{F}_s\)と\(\mathbf{G}_p\)のバイプロットなら $$ \mathbf{p} = \mathbf{D}^{-1} \mathbf{G’}_s \mathbf{x} $$ \(\mathbf{F}_p\)と\(\mathbf{G}_s\)のバイプロットなら $$ \mathbf{p} = \mathbf{G’}_s \mathbf{x} $$ となる。
4. 対応分析
[略]
5. 環境科学的な解釈
[略]
6. 応用
[略]
7. 結論
というわけで、主成分分析・対応分析のバイプロットにおける補助変数の最適な表示の仕方について考えてまいりました。
解釈のため、補助変数と主成分軸の相関をわざわざ計算する人が多いけど、実際にはそれは(適切なバイプロットを描いていれば)バイプロット上で表現されているわけだ。
云々。
—–
この論文の主眼はどっちかというと4節、対応分析における補助変数・補助点の表示にあると思うんだけど(そしてそのなかにはなんらかオリジナリティがあるんだろうけど)、いまちょっと忙しいもんで(嘘です。基本的に日々ぼーっとしてます)、申し訳ないけど読まずに飛ばした。いずれ必要になったら読みましょー。
まあとにかく、あれだ、主成分分析の場合、補助変数のベクトルは補助変数の主成分への回帰係数を表しているわけだ。はいはい、了解しましたです。