主成分の求め方　　　　　Last modified: Aug 26, 2015

  1. $p$ 個の変数を $X_{1}, X_{2}, \dots , X_{p}$，これらの重み付け合成変量を $Z_{1}, Z_{2}, \dots , Z_{m}$ とする（ $m \leqq p$ ）。

    \[ \left \{ \begin{align*} Z_1 &= L_{11}\ X_1 + L_{12}\ X_2 + \dots + L_{1p}\ X_p \\ &\vdots \\ Z_i &= L_{i1}\ X_1 + L_{i2}\ X_2 + \dots + L_{ip}\ X_p \\ &\vdots \\ Z_m &= L_{m1}\ X_1 + L_{m2}\ X_2 + \dots + L_{mp}\ X_p \\ \end{align*} \right . \] ただし，$L_{i1}^2 + L_{i2}^2 + \dots + L_{ip}^2 = 1,\ (i=1, 2, \dots, m)$ とする。

  2. このような $m$ 個の合成変量において，以下のような性質を持つものを考える。

  3. $p$ 個の変数の相関係数行列の固有値を $\lambda_{1} \geqq \lambda_{2} \geqq \dots \geqq \lambda_{i} \geqq \dots \geqq \lambda_{m} \geqq \dots \geqq \lambda_{p} \geqq 0$ としたとき，$\lambda_{i}$ に対応する固有ベクトルを重みとした合成変量が $Z_{i}$ に対応し，$Z_{i}$ の分散が $\lambda_{i}$ に等しくなる（ 固有ベクトルは互に直交する，すなわち，互に相関が 0 である ）。

    > d <- iris[,1:4]
> head(d)
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4
> p <- ncol(d) # 変数の個数

> ( r <- cor(d) ) # 相関係数
             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000

> ans <- eigen(r) # 固有値・固有ベクトルを求める

> ( lambda <- ans$values ) # 固有値
[1] 2.91849782 0.91403047 0.14675688 0.02071484

> sum(lambda) # 固有値の和は変数の個数に等しいことを確認
[1] 4

> ( L <- ans$vectors ) # 固有ベクトル（列方向）
           [,1]        [,2]       [,3]       [,4]
[1,]  0.5210659 -0.37741762  0.7195664  0.2612863
[2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
[3,]  0.5804131 -0.02449161 -0.1421264 -0.8014492
[4,]  0.5648565 -0.06694199 -0.6342727  0.5235971

> round(t(L)%*%L, 10) # 固有ベクトルが直交することを確認
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1
  4. $Z_{1}, Z_{2}, \dots , Z_{m}$ は主成分と呼ばれ，そのうちで最も分散の大きい $Z_{1}$ は第 1 主成分，次に分散の大きい $Z_{2}$ は第 2 主成分，以下順に第 $m$ 主成分と呼ばれる。

    > d2 <- sapply(d, scale) # 標準化したデータ
> Z <- d2 %*% L # 標準化したデータに重み L を掛けて主成分（主成分得点）を得る
> head(Z)
          [,1]       [,2]        [,3]         [,4]
[1,] -2.257141 -0.4784238  0.12727962  0.024087508
[2,] -2.074013  0.6718827  0.23382552  0.102662845
[3,] -2.356335  0.3407664 -0.04405390  0.028282305
[4,] -2.291707  0.5953999 -0.09098530 -0.065735340
[5,] -2.381863 -0.6446757 -0.01568565 -0.035802870
[6,] -2.068701 -1.4842053 -0.02687825  0.006586116
  :
> apply(Z, 2, var) # 各列（主成分）の分散が固有値に等しいことの確認
[1] 2.91849782 0.91403047 0.14675688 0.02071484

> head(ans2$x) # prcomp の結果の x 要素と上述の Z が等しいことの確認
           PC1        PC2         PC3          PC4
[1,] -2.257141 -0.4784238  0.12727962  0.024087508
[2,] -2.074013  0.6718827  0.23382552  0.102662845
[3,] -2.356335  0.3407664 -0.04405390  0.028282305
[4,] -2.291707  0.5953999 -0.09098530 -0.065735340
[5,] -2.381863 -0.6446757 -0.01568565 -0.035802870
[6,] -2.068701 -1.4842053 -0.02687825  0.006586116
  :
> all.equal(c(ans2$x), c(Z))
[1] TRUE
  5. 各主成分と，もとの各変数の間の相関係数は主成分負荷量と呼ばれる。主成分負荷量は，第 $i$ 主成分の重み $L_{i1}, L_{i2}, \dots , L_{ip}$ に，対応する固有値の平方根をかけたものである。すなわち，$Z_{i}$ と変数 $X_{1}$ の相関係数は $a_{1i} = L_{i1} \sqrt{\lambda_{i}}$，$Z_{i}$ と変数 $X_{2}$ の相関係数は $a_{2i} = L_{i2} \sqrt{\lambda_{i}}$ などとなる。

    > ( a <- t(t(L)*sqrt(lambda)) ) # 主成分負荷量
           [,1]        [,2]        [,3]        [,4]
[1,]  0.8901688 -0.36082989  0.27565767  0.03760602
[2,] -0.4601427 -0.88271627 -0.09361987 -0.01777631
[3,]  0.9915552 -0.02341519 -0.05444699 -0.11534978
[4,]  0.9649790 -0.06399985 -0.24298265  0.07535950

> a %*% t(a) # p 個の主成分負荷量の積は，相関係数行列に等しくなることを確認
           [,1]       [,2]       [,3]       [,4]
[1,]  1.0000000 -0.1175698  0.8717538  0.8179411
[2,] -0.1175698  1.0000000 -0.4284401 -0.3661259
[3,]  0.8717538 -0.4284401  1.0000000  0.9628654
[4,]  0.8179411 -0.3661259  0.9628654  1.0000000

> rowSums(a^2) # 主成分負荷量の行方向の二乗和は後述の寄与率（全ての主成分を抽出すれば 1）
[1] 1 1 1 1

> m <- 2 # 作成する主成分の個数（1 より大きい固有値の個数とするのが普通）
> ( a <- a[,1:m] ) # m 個の主成分を取り出す（m < p）
           [,1]        [,2]
[1,]  0.8901688 -0.36082989
[2,] -0.4601427 -0.88271627
[3,]  0.9915552 -0.02341519
[4,]  0.9649790 -0.06399985

> a %*% t(a) # m 個の主成分についての主成分負荷量の積は，相関係数行列の近似値になることを確認
            [,1]        [,2]       [,3]       [,4]
[1,]  0.92259864 -0.09109425  0.8911004  0.8820872
[2,] -0.09109425  0.99091932 -0.4355879 -0.3875343
[3,]  0.89110035 -0.43558792  0.9837300  0.9583285
[4,]  0.88208719 -0.38753432  0.9583285  0.9352804

> rowSums(a^2) # m 個の主成分負荷量の行方向の二乗和は後述の寄与率（1 より小さい）
[1] 0.9225986 0.9909193 0.9837300 0.9352804
  6. 主成分分析の結果を表 1 のように表す。

    表 1．主成分分析の結果の表現

    		　第1主成分　 　第2主成分　 　$\dots$　 　第$m$主成分　 　　寄与率*　　
    $X_{1}$ $a_{11}$ $a_{12}$ $\dots$ $a_{1m}$ $\sum_{k=1}^m a_{1k}^{2}$
    $X_{2}$ $a_{21}$ $a_{22}$ $\dots$ $a_{2m}$ $\sum_{k=1}^m a_{2k}^{2}$
    ： $\vdots$ $\vdots$ $\dots$ $\vdots$ $\vdots$
    $X_{p}$ $a_{p1}$ $a_{p2}$ $\dots$ $a_{pm}$ $\sum_{k=1}^m a_{pk}^{2}$
    固有値 $\sum_{j=1}^p a_{j1}^{2}$ $\sum_{j=1}^p a_{j2}^{2}$ $\dots$ $\sum_{j=1}^p a_{jm}^{2}$
    寄与率** $\displaystyle \frac{\sum_{j=1}^p a_{j1}^{2}}{p}$ $\displaystyle \frac{\sum_{j=1}^p a_{j2}^{2}}{p}$ $$\dots$$ $\displaystyle \frac{\sum_{j=1}^p a_{jm}^{2}}{p}$ 

                   第1主成分  第2主成分  寄与率
Sepal.Length       0.890     -0.361   0.923
Sepal.Width       -0.460     -0.883   0.991
Petal.Length       0.992     -0.023   0.984
Petal.Width        0.965     -0.064   0.935
固有値             2.918      0.914
寄与率            72.962%    22.851%
累積寄与率        72.962%    95.813%
  7. 寄与率* の欄は，各変数が $m$ 個の主成分でどれくらい説明されるかを表す（ $0 \leqq$ 寄与率 $\leqq 1$ ）。

    > ( Eigenvalues <- colSums(a^2) ) # 主成分負荷量の縦方向の二乗和は固有値に等しいことの確認
[1] 2.9184978 0.9140305
  8. もう一つの寄与率** はパーセントで表され，各主成分がもとの情報をどれくらい説明しているかを表すもので，相関係数から出発した主成分分析の場合には，$p$ 個の変数の持つ情報量の合計は $p$ であるので，例えば第 1 主成分の寄与率は $\displaystyle \frac{\displaystyle \sum_{j=1}^p a_{j1}^{2}}{p}$ である。 
    > Eigenvalues/p*100 # 各主成分がどれほどの情報を持っているか
[1] 72.96245 22.85076
> cumsum(Eigenvalues/p)*100 # 累積寄与率
[1] 72.96245 95.81321

演習問題：

応用問題：

・ 次のページへ進む　　 ・ 主成分分析の最初のページへ戻る
・ E-mail to Shigenobu AOKI