はじめに
千葉大学・株式会社Nospareの川久保です．今回は多変量解析の代表的な分析手法の1つである主成分分析について解説します．
主成分分析は，高次元データをできるだけ情報を損失しない形で次元を落とすことを目的にしています．応用例をいくつか簡単に紹介します．まずは画像データの圧縮です．ピクセルごとにグレースケールが観測されているデータ（ピクセル数がデータの次元）を，主成分分析により次元圧縮し，主成分ベクトルを用いて復元することが可能です．次に主成分回帰と呼ばれる分析手法があります．これは，主成分分析により次元削減された説明変数を用いて，元データよりも少ない数の説明変数を用いて回帰分析を行うというものです．またマーケティングへの応用も考えられ，ブランドや商品の特徴を，より少ない要素で説明し解釈可能性を上げる効果が期待されます．
次元削減の方法のイメージを簡単に説明します．$p$次元データ$(x_{i1},\dots,x_{ip})$が$n$個観測されているとします（$i=1,\dots,n$）．この$p$次元データに何らかの変換を施して，$m(<p)$次元データ$y_i = (y_{i1},\dots,y_{im})$を得たいとします．その何らかの変換ですが，以下のような元データ$x$の線形結合を考えます．
\begin{equation}
\begin{split}
y_{i1} &= w_{11}x_{i1} + w_{12}x_{i2} + \dots + w_{1p}x_{ip}, \\
y_{i2} &= w_{21}x_{i1} + x_{22}x_{i2} + \dots + w_{2p}x_{ip}, \\
\vdots \\
y_{im} &= w_{m1}x_{i1} + w_{m2}x_{i2} + \dots + w_{mp}x_{ip}.
\end{split}
\tag{1}
\end{equation}
ウェイト$w$をどのように決定するか，以下で解説していきます．
主成分の計算
第一主成分
まずは第一主成分，すなわち(1)式の1行目におけるウェイトのベクトルを求めます．以下のようなデータの変換を考えます．
\begin{align}
y_i &= w_1x_{i1} + w_2x_{i2} + \dots + w_px_{ip} \\
&= w^\top x_i.
\end{align}
ただし，$w = (w_1,\dots,w_p)^\top, x_i = (x_{i1},\dots,x_{ip})^\top$で，$\| w \|^2 = w^\top w = 1$とします．このとき$y_i$は，$x_i$の$w$への正射影ベクトルの符号付き長さ（影の長さ）とみなせます．つまり以下の図における$x_i$から射影軸$w$におろした垂線の足（$H$）と原点$O$との距離で，射影軸の目盛を読むイメージです．ここで射影軸上で最もデータがばらつくように（つまり$y_1,\dots,y_n$の分散が最大になるように）$w$を決定し，その$w$のことを第一主成分と呼びます．射影については，私の過去の記事を参照してください．
射影軸上のデータ$y_1,\dots,y_n$の分散は，
\begin{align}
s_y^2 &= \frac{1}{n} \sum_{i=1}^n (y_i - \bar{y})^2 \\
&= \frac{1}{n}\sum_{i=1}^n \{ w^\top (x_i - \bar{x}) \}^2 \\
&= w^\top \left\{ \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(x_i - \bar{x})^\top \right\} w \\
&= w^\top S w.
\end{align}
ただし$S$は元の$p$次元データ$x$の標本共分散行列です．ゆえに，分散$s_y^2$が最大となる射影軸$w$を求める問題は，$w^\top w = 1$の制約のもとでの$w^\top S w$の最大化問題に帰着されます．最大化問題の解は，$S$の最大固有値に対応する固有ベクトルに他ならないのですが，そのことを以下で示していきます．
共分散行列$S$の固有値$\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_p \geq 0$，それに対応する固有ベクトルを$v_1,v_2,\dots,v_p$とすると，$V = [v_1,\dots,v_p]$と$\Lambda = \mathrm{diag}(\lambda_1,\dots,\lambda_p)$を用いて，
$$
S = V \Lambda V^\top
$$
と書けます．ここで$V$は$p \times p$の直交行列，すなわち$V^\top V = VV^\top = I_p$が成り立ちます．ここで，$\| w \| = 1$を満たす$w^\top S w$の最大化問題について，
\begin{equation}
\begin{split}
w^\top S w &= w^\top V \Lambda V^\top w \\
&= \tilde{w}^\top \Lambda \tilde{w} \quad (\tilde{w} := V^\top w) \\
&= \sum_{j=1}^p \tilde{w}_j^2 \lambda_j
\end{split}
\tag{2}
\end{equation}
なる式変形が成り立ちますが，$\| w \| = 1$のもとでは，$\| \tilde{w} \|^2 = w^\top VV^\top w = w^\top w = 1$となります．よって，今考えている最大化問題は，(2)式の最右辺を$\tilde{w}_1^2 + \dots + \tilde{w}_p^2 = 1$の制約のもとで最大化する問題と同値になり，その最大値は$\lambda_1$が達成されることが容易にわかります．つまり第一主成分による変換データの分散は$\lambda_1$です．また最大化元$w$は，固有値$\lambda_1$に対応する固有ベクトル$v_1$であり，これが求めていた分散を最大にする射影ベクトルです．つまり元の$p$次元データの標本共分散行列$S$の最大固有値に対応する固有ベクトル$v_1$が，第一主成分であることがわかりました．
第二主成分以降の導出
実は第二主成分以降も，それぞれ第$j$固有値に対応する固有ベクトルが第$j$主成分となります．第二主成分は，第一主成分$v_1$に直交（$w ^\top v_1 = 0$）しているという条件のもとで，$w^\top S w$（変換データの分散）を最大とする$\| w \| = 1$を満たすベクトル$w$です．このとき，
\begin{equation}
\begin{split}
w^\top S w &= w^\top V \Lambda V^\top w \\
&= w^\top [v_1, v_2,\dots,v_p]
\begin{bmatrix}
\lambda_1 & & & \\
& \lambda_2 & & \\
& & \ddots & \\
& & & \lambda_p
\end{bmatrix}
\begin{bmatrix}
v_1^\top \\
v_2^\top \\
\vdots \\
v_p^\top
\end{bmatrix}
w \\
&= [0, w^\top v_2,\dots,w^\top v_p]
\begin{bmatrix}
\lambda_1 & & & \\
& \lambda_2 & & \\
& & \ddots & \\
& & & \lambda_p
\end{bmatrix}
\begin{bmatrix}
0 \\
v_2^\top w \\
\vdots \\
v_p^\top w
\end{bmatrix}
\quad (\because v_1^\top w = 0) \\
&= \sum_{j=2}^p \tilde{w}_j^2 \lambda_j \quad (\tilde{w} = (\tilde{w}_1,\dots,\tilde{w}_p)^\top := V^\top w)
\end{split}
\tag{3}
\end{equation}
なる式変形が成り立ち，先ほどの議論と同様，$\| w \| = 1$のもとでは$\| \tilde{w} \| = 1$が成り立ちます．よって，今考えている最大化問題は，(3)式の最右辺を$\tilde{w}_1^2 + \dots + \tilde{w}_p^2 = 1$のもとで最大化する問題と同値で，その最大値は$\lambda_2$が達成されることがわかります．またその最大化元$w$は固有値$\lambda_2$に対応する固有ベクトル$v_2$であり，これが第二主成分です．
第三主成分以降は，それまでの主成分（第三主成分を求める問題の場合は，第一主成分と第二主成分）とすべて直交し（$w^\top v_1 = w^\top v_2 = 0$），かつ$\| w \| = 1$のもとで$w^\top S w$を最大とするベクトル$w$として求めます．(3)式と同様の式変形を行うことで，これが第$j$固有値に対応する固有ベクトルであることは容易に確認できます．
2次元データにおける主成分分析のイメージは以下の図のようになります．
主成分得点と寄与率
第$j$主成分$w_j = (w_{j1},w_{j2},\dots,w_{jp})^\top (= v_j)$を固有値問題により求めることができれば，$i$番目のデータを変換した以下の値
$$
y_{ij} = w_{j1}x_{i1} + w_{j2}x_{i2} + \dots + w_{jp}x_{ip} = w_j^\top x_i
$$
を計算することができます．この値のことを，$i$番目のデータの第$j$主成分得点と呼びます．主成分得点の大小から，それぞれのデータの第$j$主成分における傾向や意味を見出すのに役立ちます．
また第$j$主成分で変換したデータ$y_{1j},\dots,y_{nj}$の分散が$\lambda_j$であったことに留意すると，以下の
$$
\frac{\lambda_j}{\lambda_1 + \lambda_2 + \dots + \lambda_p}
$$
の値は，全体のうち第$j$主成分が持っている情報の割合と理解できます．これを寄与率と呼びます．また，
$$
\frac{\lambda_1 + \lambda_2 + \dots + \lambda_j}{\lambda_1 + \lambda_2 + \dots + \lambda_j + \dots + \lambda_p}
$$
は，第$j$主成分までの累積寄与率と呼ばれます．逆に，1から累積寄与率を引いた値が，第$j$主成分までで次元削減をした場合の，情報損失の割合と理解できます．
数値例
簡単な数値例を見てみます．5教科（英・数・国・理・社）の5次元成績データを，主成分分析により次元削減し，それぞれの主成分がどのような能力と解釈できるかを検討します．
5教科の得点を表す変数をそれぞれ$x_1,x_2,x_3,x_4,x_5$とし，標本共分散行列が以下のように求まったとします．
S =
\begin{bmatrix}
220 & 12 & 156 & 10 & 222 \\
12 & 286 & 63 & 240 & 16 \\
156 & 63 & 164 & 31 & 153 \\
10 & 240 & 31 & 250 & 17 \\
222 & 16 & 153 & 17 & 248
\end{bmatrix}
この$S$の固有値・固有ベクトルを計算します．Rではeigenという関数が実装されています．
> PCA <- eigen(S)
> PCA
eigen() decomposition
$values
[1] 606.584273 479.946697 56.047456 16.848549 8.573024
$vectors
[,1] [,2] [,3] [,4] [,5]
[1,] 0.5063307 -0.3475368 -0.1391771 -0.1777353 0.7562323
[2,] 0.3849703 0.6210881 0.2989299 0.5738915 0.2175703
[3,] 0.4419006 -0.1614720 0.7282486 -0.3778524 -0.3248572
[4,] 0.3372365 0.5852321 -0.4436130 -0.5636948 -0.1709696
[5,] 0.5351845 -0.3534089 -0.4051323 0.4225333 -0.4959970
PCA$valuesが$S$の5つの固有値で，PCA$vectorsの各列が$S$の固有ベクトルたちです．第二主成分までの累積寄与率を計算してみます．
> la <- PCA[[1]]
> sum(la[1:2]) / sum(la)
[1] 0.9302491
第二主成分までで，データ全体のおよそ93パーセントの情報を持っていると理解できます．また固有ベクトル（主成分）に着目すると，第一主成分ベクトルはすべて正の値であることから，これは総合的な能力と理解することができます．また，第二主成分は，$x_2$と$x_4$（数学と理科）に対応するウェイトが正の値，残りが負の値であることから，理数系の能力を表すと解釈することができそうです．
おわりに
株式会社Nospareには，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．
