主成分分析と言うデータ分析手法がある。複数の評価項目のあるデータの全体像をつかむ際に、その評価項目の特性を分析して、より小数の新たな評価項目を見出し、データの全体像をより簡単に把握、可視化しようという試みだ。
例えば、中学生の成績を評価するとしよう。いわゆる5大教科（国、数、英、理、社）に加えて体育、音楽、技術など、それぞれの評価項目から全体から把握するよりも、これらに起因するより簡単な評価項目（評価軸）を見出すことができれば、おぼろげながらも、より簡単な把握が可能になる。例えば、高校入試科目とそれ以外とか、筆記試験を伴う科目と実技試験を伴う科目とか。このような評価軸を統計的に見出だし、リサーチャーが解釈し、データの傾向を把握するのが主成分分析だ。
ここでは主成分分析そのものを紹介するのではなく、Rを用いてその分析手順を紹介していく。分析対象データはRに含まれているものを用いるので、RStudioだけで手順を再現することができる。ここでは次の環境を用いている。
- R 3.5.2
- RStudio 1.2.5001
ちなみに現時点でのR、最新版は3.6.2だ。RStudioと共に参照先は、この投稿末尾にまとめている。
ライブラリとデータの準備
この投稿では、次のライブラリを用いる。
| FactoMineR | 主成分分析に、関数PCAを用いる。 |
| Factoextra | 主成分分析の視覚化に用いる。 |
| corrplot | 相関行列の視覚化に用いる。 |
| LifeCycleSavings | 今回の分析対象データ。 |
いずれについても、この投稿末尾に参照先をまとめている。
主成分分析に用いる関数はPCA以外にもprcompやprincompなどを用いる場合もあるが、ここではPCAを用いている。端的には、私が気に入っているから、と言うのがその理由だ。特に分析結果を構造的に出力してくれる点を気に入っている。例えば「mypca = PCA(...)」の出力結果の一部は、次のような具合だ。
| mypca$eig | 固有値 | ||
| mypca$var | 変数（評価項目) | $coord | 座標 |
| $cor | 主成分との相関 | ||
| $cos2 | スコア | ||
| $contrib | 寄与率 | ||
| mypca$ind | 個々のデータ（分析対象） | $coord | 座標 |
| $cos2 | スコア | ||
| $contrib | 寄与率 | ||
| mypca$$ind.sup | 個々のデータ（予測対象） | $coord | 座標 |
| cos2 | スコア |
LifeCycleSavingsは可処分所得に対する貯蓄率のデータだ。国別に集計されている。変数は次のように定義されている。数値は、1960～1970年の平均だ。
| sr | 個人の貯蓄率 |
| pop15 | 15歳未満の人口割合 |
| pop75 | 75歳以上の人口割合 |
| dpi | 可処分所得に対する金額 |
| ddpi | dpiの成長率 |
データの読み込みを含め、これらを利用するためのコードは次のようになる。
R code
library(FactoMineR) library(factoextra) library(corrplot) mydata = LifeCycleSavings
主成分の解釈
まず予測無しに、主成分分析の流れを追ってみる。
関数PCAを用いて、mydataに対する主成分分析を実施する。分析結果をmypcaに格納する。
mypca$eigで、分析結果から各主成分の固有値を表示させると、
R code
mypca = PCA(mydata, scale.unit = TRUE, graph = FALSE) mypca$eig
Outputs
eigenvalue percentage of variance cumulative percentage of variance comp 1 2.82207781 56.441556 56.44156 comp 2 1.25606655 25.121331 81.56289 comp 3 0.60452546 12.090509 93.65340 comp 4 0.23964496 4.792899 98.44630 comp 5 0.07768522 1.553704 100.00000
「cumulative percentage of variance」を参照すると、第１、２主成分で、データ全体の81.6%を説明できることが分かる。
fviz_pca_varで各変数と主成分との相関を確認する。スコア（cos2）、寄与率（contrib）と共に確認してみる。矢印の長さが主成分との相関の表している。矢印が長ければ（円に近ければ）、相関は高い。pop15～dpiは第１主成分との相関が高く、ddpiが第２主成分との相関が高いことが分かる。
ちなみに第１主成分が横軸、第２主成分が縦軸だ。
R code
fviz_pca_var(mypca, col.var = "cos2", gradient.cols = c("blue", "red"), repel = TRUE) fviz_pca_var(mypca, col.var = "contrib", gradient.cols = c("blue", "red"), repel = TRUE)
具体的な数値で把握するため、corrplotを用いると、第１、２主成分が形成する４象限を解釈しやすい。
R code
corrplot(mypca$var$cos2, addCoef.col = "gray") corrplot(mypca$var$contrib, is.corr = FALSE, addCoef.col = "gray")
第１主成分左側はpop15、右側はpop75、dpiの領域だ。第1主成分の左から右にかけて可処分所得が大きくなる。左へ行くほど若い人の領域、右に行くほど年寄りの領域、と解釈できる。
第２主成分上側はddpiの領域だ、上に行くほど可処分所得の成長率が高く、下に行くほど低いと解釈できる。
主成分の解釈はリサーチャーに委ねられる。今回のデータは、たまたま素直な出力を返してくれたので、主成分の解釈はかなり率直だった。データによっては解釈に苦しむ出力を返すこともある。考えることが好きな人には面白いところでもある。
個々のデータのプロット
まず主成分上のすべてのデータをプロットしてみる。fviz_pca_indを用いる。
R code
fviz_pca_ind(mypca, repel = TRUE)
大変、見難い。各データを第１、２主成分に反映するにあたり、元のデータから失われている情報がある。各データのスコアを参照し、あまり低いデータは表示しないよう、フィルタをかける。まずスコアを確認する。
R code
fviz_cos2(mypca, choice = "ind", axes = 1:2)
第１、２象限に対するcos2について、0.5を超えるあたりで一区切りするのが妥当に見えるが、それでもまだデータ数は多い。ここではcos2 > 0.8を表示対象にしてみる。かなりすっきりした。
R code
fviz_pca_ind(mypca, repel = TRUE, select.ind = list(cos2 = 0.8))
全体の解釈
fviz_pca_biplotで変数と、個々のデータをプロットしてみる。注意しなければならないのは、個々のデータのポジションと、変数を示す矢印の方向に注目することだ。変数を示す矢印と同じ方向にあるポジションが、変数の傾向を色濃く反映しているデータであり、反対方向にあれば傾向は弱くなる。
R code
fviz_pca_biplot( mypca, col.var = "blue", col.ind = "cos2", gradient.cols = c("blue", "red"), repel = TRUE, select.ind = list(cos2 = 0.8))
例えばChinaは若い領域におり、可処分所得の金額は控えめなものの、その成長率は高いと解釈できる。60～70年代、まだ中国は今のような超大国ではなかった。
Japanは年寄りの領域におり、可処分所得は高く、その成長率も高い。実際、60～70年代と言えば高度経済成長期だった。
可処分所得が高いものの、成長率の低い領域には北欧、西欧諸国が確認できる。一方、南米諸国は可処分所得が低く、成長率も低いことが分かる。
予測
貯蓄率と日本を除外したデータで主成分分析し、その結果をもとに、それらがどのようにプロットされるのかを確認してみる。つまり、主成分分析の結果に基づいて、それらを予測するということだ。
貯蓄率はデータの1列目、日本はデータの23行目に記録されている。主成分分析の再実施に際し、日本をind.supに、貯蓄率をquanti.supに指定する。
R code
mypca = PCA(mydata, ind.sup = 23, quanti.sup = 1 ,scale.unit = TRUE, graph = FALSE)
PCA実行後、ind.supとquanti.supを確認するとともに、それらをfviz_pca_indにプロットしてみる。
R code
mypca$ind.sup mypca$quanti.sup
Outputs
$coord Dim.1 Dim.2 Dim.3 Dim.4 Japan 0.3953987 1.620821 0.09006812 -0.8875493 $cos2 Dim.1 Dim.2 Dim.3 Dim.4 Japan 0.04367951 0.7339682 0.002266467 0.2200859 $dist Japan 1.891892 $coord Dim.1 Dim.2 Dim.3 Dim.4 sr 0.3629719 0.275223 -0.1636127 -0.2067154 $cor Dim.1 Dim.2 Dim.3 Dim.4 sr 0.3629719 0.275223 -0.1636127 -0.2067154 $cos2 Dim.1 Dim.2 Dim.3 Dim.4 sr 0.1317486 0.07574771 0.0267691 0.04273127
R code
myplot = fviz_pca_biplot( mypca, col.var = "orange", col.ind = "cos2", gradient.cols = c("gray", "green"), repel = TRUE, select.ind = list(cos2 = 0.8)) myplot = fviz_add(myplot, mypca$ind.sup$cos2, color = "blue", repel = TRUE) myplot
分析結果であるmypcaのindには予測対象が含まれないため、fviz_addでind.supのプロットを指定し、myplotと重ね書きする必要がある。
先ほどとは異なる位置に、貯蓄率srが点線の青矢印で表示され、日本が青でプロットされているのが分かる。
参照
www.r-project.org
rstudio.com
www.rdocumentation.org
www.rdocumentation.org
www.rdocumentation.org
www.rdocumentation.org
www.rdocumentation.org
コード全文
